NASA— CR- 19956 9 


/=> // /) (~ 


3 \ iT7C 

Sensitivity Analysis for Aeroacoustic and Aeroelastic 

Design of Turbomachinery Blades j)^> 

Christopher B. Lorence 
Graduate Research Assistant 


and 


Kenneth C. Hall 
Principal Investigator 

Department of Mechanical Engineering and Materials Science 

School of Engineering 
Duke University 
Durham, NC 27708-0300 



DUKE 


SCHOOL OF 
ENGINEERING 


A Final Technical Report on Research Conducted Under 
NASA Grant NAG3-1433 

(Grant Title: Aeroacoustic Sensitivity Analysis and Optimal Aeroacoustic 

Design of Turbomachinery Blades) 


Submitted to NASA Lewis Research Center 


/« uuL'iiuuvu lv/ nx lux i juv vy 

c May 1995 


- i (NIP 3 ^95-0 0 590)" SENSITIVITY 
ANALYSIS FOR AEROACOUSTIC AND 
AEROELASTIC DESIGN OF 
TURBOMACHINERY BLADES Final Report 
(Cuke Unlv.) 162 p 


N96-13226 


Unci as 


G3/07 0073244 



Abstract 


A new method for computing the effect that small changes in the airfoil shape and 
cascade geometry have on the aeroacoustic and aeroelastic behavior of turbomachin- 
ery cascades is presented. The nonlinear unsteady flow is assumed to be composed of 
a nonlinear steady flow plus a small perturbation unsteady flow that is harmonic in 
time. First, the full potential equation is used to describe the behavior of the non- 
linear mean (steady) flow through a two-dimensional cascade. The small disturbance 
unsteady flow through the cascade is described by the linearized Euler equations. Us- 
ing rapid distortion theory, the unsteady velocity is split into a rotational part that 
contains the vorticity and an irrotational part described by a scalar potential. The 
unsteady vorticity transport is described analytically in terms of the drift and stream 
functions computed from the steady flow. Hence, the solution of the linearized Euler 
equations may be reduced to a single inhomogeneous equation for the unsteady poten- 
tial. The steady flow and small disturbance unsteady flow equations are discretized 
using bilinear quadrilateral isoparametric finite elements. The nonlinear mean flow 
solution and streamline computational grid are computed simultaneously using New- 
ton iteration. At each step of the Newton iteration, LU decomposition is used to solve 
the resulting set of linear equations. The unsteady flow problem is linear, and is also 
solved using LU decomposition. Next, a sensitivity analysis is performed to deter- 
mine the effect small changes in cascade and airfoil geometry have on the mean and 
unsteady flow fields. The sensitivity analysis makes use of the nominal steady and 
unsteady flow LU decompositions so that no additional matrices need to be factored. 
Hence, the present method is computationally very efficient. To demonstrate how the 
sensitivity analysis may be used to redesign cascades, a compressor is redesigned for 
improved aeroelastic stability and two different fan exit guide vanes are redesigned 
for reduced downstream radiated noise. In addition, a framework detailing how the 
two-dimensional version of the method may be used to redesign three-dimensional 
geometries is presented. 


Contents 


1 Introduction 8 

1.1 Statement of the Problem 8 

1.1.1 Aeroelastic Problems in Turbomachines 9 

1.1.2 Aeroacoustic Problems in Turbomachines 11 

1.2 Previous Work 12 

1.2.1 Unsteady Aerodynamic Models 12 

1.2.2 Design Methods and Optimization 15 

1.3 Present Method 16 

1.4 Outline of Report 17 

2 Nominal Flow Field Description 18 

2.1 Euler Equations 18 

2.1.1 Goldstein Decomposition 20 

2.1.2 Atassi Decomposition 24 

2.2 Extension of Bateman’s Variational Principle 27 

2.2.1 Nonlinear Full Variational Principle 28 

2.2.2 Steady Flow Variational Principle 29 

2.2.3 Unsteady Flow Variational Principle With Deforming Grid and 

Vorticity 29 

2.3 Near-Field Boundary Conditions 34 

2.3.1 Periodic Boundary Condition 34 

2.3.2 Airfoil Surface Boundary Condition 36 

2.3.3 Wake Boundary Condition 37 

2.4 Far-Field Boundary Conditions 38 

2.4.1 Steady Flow 38 

2.4.2 Unsteady Flow 39 

3 Numerical Solution Method 44 

3.1 Steady Flow Solution Procedure 44 

3.1.1 Grid Generation 44 

3.1.2 Airfoil Definition and Spline Notation 48 

3.1.3 Finite Element Discretization 51 

3.1.4 Near-Field Boundary Conditions 53 

3.1.5 Far-Field Boundary Conditions 57 

3.1.6 Assembly and Solution 58 


1 


2 


CONTENTS 


3.2 Unsteady Flow Solution Procedure 61 

3.2.1 Drift Function and Rotational Velocity Calculation 61 

3.2.2 Unsteady Grid Generation 63 

3.2.3 Finite Element Discretization 64 

3.2.4 Near-Field Boundary Conditions 66 

3.2.5 Far-Field Boundary Conditions 70 

3.2.6 Assembly and Solution 75 

4 Sensitivity Analysis 76 

4.1 General Procedure 76 

4.2 Steady Flow Sensitivity Analysis 79 

4.2.1 Prescribing the Perturbation 79 

4.2.2 Near-Field Boundary Conditions 80 

4.2.3 Far-Field Boundary Conditions 82 

4.2.4 Sensitivity of the Stream Function 83 

4.3 Unsteady Flow Sensitivity Analysis 84 

4.3.1 Sensitivity of the Drift Function 84 

4.3.2 Sensitivity of the Rotational Velocity 84 

4.3.3 Sensitivity of the Blade Structural Dynamics 85 

4.3.4 Sensitivity of the Grid Motion 86 

4.3.5 Sensitivity of the Finite Element Assembly 86 

4.3.6 Near-Field Boundary Conditions 87 

4.3.7 Far-Field Boundary Conditions 88 

5 Results 89 

5.1 Aeroelastic Analysis and Design of a Compressor 89 

5.1.1 Steady Flow Through a Compressor 89 

5.1.2 Unsteady Flow Through a Compressor 90 

5.1.3 Sensitivity Analysis 95 

5.1.4 Redesign of a Compressor for Aeroelastic Stability 100 

5.1.5 Computational Efficiency 105 

5.2 Aeroacoustic Analysis and Design of a Fan Exit Guide Vane 107 

5.2.1 Steady Flow Through a Fan Exit Guide Vane 108 

5.2.2 Unsteady Flow Through a Fan Exit Guide Vane 109 

5.2.3 Sensitivity Analysis 110 

5.2.4 Redesign of an EGV for Reduced Acoustic Response 115 

5.2.5 Computational Efficiency 125 

5.3 Modern Fan Exit Guide Vane 125 

5.3.1 Nominal Analysis 126 

5.3.2 Sensitivity Analysis 127 

5.3.3 Redesign of EGV for Reduced Acoustic Response 129 

5.4 Summary 132 


CONTENTS 


3 


6 Application to Three-Dimensional Problems 134 

6.1 Acoustic Modes in an Annular Duct 134 

6.2 Calculation of Far-Field Unsteady Pressure 138 

6.3 Calculation of the Outgoing Pressure 139 

7 Conclusions and Future Work 143 

7.1 Conclusions 143 

7.2 Future Work 144 

7.2.1 Other Flow Models 144 

7.2.2 Multidisciplinary Optimization 145 

7.2.3 Multiple Blade Rows 145 

7.2.4 Three-Dimensional Problems 146 

A Nomenclature 148 

B Sensitivity of the NACA Modified Four-Digit Airfoil Definition 152 


Bibliography 


155 



List of Figures 


1.1 Axial compressor or fan characteristic map showing principal types of 

flutter and regions of occurrence (adapted from [3]) 10 

1.2 Typical compressor resonance (or Campbell) diagram (adapted from 

[5]) 11 

2.1 Contours of drift function, A (x,y), and stream function, T(x,y), for a 

typical fan exit guide vane 22 

2.2 Left, contours of the stream function, T, near a solid plane boundary. 
Right, contours of the drift function, A, near an airfoil stagnation point. 23 

2.3 Contours of computational coordinates (£, 77) for undeformed (top) and 

deformed (bottom) grids for a cascade of airfoils. Airfoils are pitching 
about their midchords with an interblade phase angle, a, of 180°. . . 31 

2.4 Locations of boundaries for a typical cascade 35 

2.5 Coordinate systems for downstream far-field analysis 42 

3.1 Typical computational grid for a fan exit guide vane 45 

3.2 Typical computational grid in (E,\?) coordinates. The computational 

node numbering convention is also illustrated 47 

3.3 Typical distribution of boundary grid points 48 

3.4 Detail of airfoil spline definition. Arrows indicate direction of airfoil 

spline 49 

3.5 Typical computational grid illustrating modifications for airfoils with 

blunt trailing edges. Multiple passages shown for clarity. The high- 
lighted lines define the actual computational grid 50 

4.1 Overview of quantities in sensitivity analysis 77 

5.1 Steady surface pressure of cascade of NACA 5506 airfoils. M_oo = 0.5, 

0_ oo = 55°, 0 = 45° 90 

5.2 Aerodynamic damping of cascade of NACA 5506 airfoils vibrating in 

plunge at frequencies of 0.4, 0.8, and 1.6 for a range of interblade phase 
angles 92 

5.3 Aerodynamic damping of cascade of NACA 5506 airfoils pitching about 

their midchords at frequencies of 0.4, 0.8, and 1.6 for a range of in- 
terblade phase angles 93 


4 


LIST OF FIGURES 


5 


5.4 Imaginary part of unsteady surface pressure of NACA 5506 airfoils 

pitching about their midchords. u> — 0.4, a = 60° 94 

5.5 Sensitivity of steady surface pressure of cascade of NACA 5506 airfoils 
to perturbations in thickness, camber, stagger, gap, and reflex. M_ 00 = 

0.5, fi_oo = 55°, 0 = 45° 96 

5.6 Real part of sensitivity of unsteady surface pressure of NACA 5506 air- 

foils pitching about their midchords due to perturbations in thickness, 
camber, stagger, gap, reflex, and frequency. u> = 0.4, a — 60° 98 

5.7 Imaginary part of sensitivity of unsteady surface pressure of NACA 
5506 airfoils pitching about their midchords due to perturbations in 


thickness, camber, stagger, gap, reflex, and frequency, u = 0.4, a — 60°. 99 
5.8 Steady surface pressure of cascade of redesigned airfoils. M - = 0.5, 

oo = 55° 102 


5.9 Real and imaginary parts of unsteady surface pressure of redesigned 
airfoils (Redesign A) pitching about their midchords. u> = 0.4, a = 60°. 103 

5.10 Real and imaginary parts of unsteady surface pressure of redesigned 
airfoils (Redesign B) pitching about their midchords, u = 0.4, cr = 60°. 104 

5.11 Accuracy of sensitivity analysis for various perturbation amplitudes. . 105 

5.12 Aerodynamic damping of cascade of redesigned airfoils (Redesign A) 

pitching about their midchords at a frequency of 0.4 for a range of 
interblade phase angles 106 

5.13 Aerodynamic damping of cascade of redesigned airfoils (Redesign B) 

pitching about their midchords at a frequency of 0.4 for a range of 
interblade phase angles 107 

5.14 Schematic showing wake/EGV interaction. Inserts in blade passage 

show contours of “drift” (top) and a streamline computational grid 
(bottom) 108 

5.15 Steady surface pressure for cascade of NACA 8508-65 airfoils. 0 = 

16°, n_oo = 30°, M_oo = 0.5, G = 1.0 109 

5.16 Unsteady surface pressure of cascade of NACA 8508-65 airfoils due to 

incoming vortical gust at lxBPF. t o = 3.7687, g = —144° Ill 

5.17 Unsteady surface pressure for cascade of NACA 8508-65 airfoils due 

to incoming vortical gust at 2xBPF. u = 7.54, a — —288° 112 

5.18 Contours of magnitude of unsteady pressure for cascade of NACA 
8508-65 airfoils due to incoming vortical gust at lxBPF. u = 3.77, 

g = -144° 113 

5.19 Contours of magnitude of unsteady pressure for cascade of NACA 
8508-65 airfoils due to incoming vortical gust at 2xBPF. u = 7.54, 

a = -288° 114 

5.20 Sensitivity of unsteady surface pressure on NACA 8508-65 airfoils due 
to perturbations in thickness, camber, stagger, gap, and reflex, u = 

7.54, a = —288°. □ , suction surface; A, pressure surface 116 

5.21 Sensitivity of unsteady surface pressure on NACA 8508-65 airfoils due 
to perturbations in thickness, camber, stagger, gap, and reflex, u = 

7.54, g = —288°. □, suction surface; A, pressure surface 


117 



6 


LIST OF FIGURES 


5.22 Steady surface pressure for cascade of redesigned airfoils. = 0.5, 

ft-co = 30° 119 

5.23 Real and imaginary parts of unsteady surface pressure for redesigned 

airfoils due to incident vortical gust at 2xBPF. ui = 7.54, a = —317°. 120 

5.24 Accuracy of sensitivity analysis for various perturbation amplitudes. . 121 

5.25 Contours of unsteady pressure for cascade of redesigned airfoils due to 

incoming vortical gust at lxBPF. u> = 3.77, a = —158.4° 122 

5.26 Contours of unsteady pressure for cascade of redesigned airfoils due to 

incoming vortical gust at 2xBPF. u = 7.54, a = —317° 123 

5.27 Contours of magnitude of unsteady pressure for radial station near the 

hub of a modern fan exit guide vane due to incoming vortical gust at 
2xBPF. M_oo = 0.32, G = 0.46, u = 12.5, cx = -288° 126 

5.28 Contours of magnitude of unsteady pressure for radial station near the 
midspan of a modern fan exit guide vane due to incoming vortical gust 

at 2xBPF. Af_oo = 0.56, G = 0.71, w = 7.7, cx = -288° 127 

5.29 Contours of magnitude of unsteady pressure for radial station near the 

tip of a modern fan exit guide vane due to incoming vortical gust at 
2xBPF. = 0.40, G = 0.90, to = 10.4, a = -288° 128 

5.30 Contours of magnitude of unsteady pressure for radial station near 
the hub of a redesigned modern fan exit guide vane due to incoming 
vortical gust at 2xBPF. M_ = 0.32, G = 0.46, u> = 12.5, a = —288°. 130 

5.31 Contours of magnitude of unsteady pressure for radial station near the 

midspan of a redesigned modern fan exit guide vane due to incoming 
vortical gust at 2xBPF. M - «, = 0.56, G = 0.71, u — 7.7, cr = —288°. 130 

5.32 Contours of magnitude of unsteady pressure for radial station near the 

tip of a redesigned modern fan exit guide vane due to incoming vortical 
gust at 2xBPF. Af_ c 0 = 0.40, G = 0.90, cu = 10.4, cx = -288° 133 

6.1 Annular duct geometry 135 

6.2 Typical radial mode shapes fi mn for annular duct. Hub-to-tip ratio, 
rH/ r T — 0-5; solid line, n = 0; dashed line, n = 1; dotted line, n = 2. 

Left, Fourier mode m = 1; center, Fourier mode m = 2; right, Fourier 
mode m = 8 137 



List of Tables 


5.1 Sensitivity of steady forces and moment. The nominal steady lift, 

L, is 0.2907, the nominal drag, D , is —0.0177, the nominal moment 
about the leading edge, Mle , is —0.1215, and the nominal lift in the 
^-direction, Ly, is 0.1931 97 

5.2 Sensitivity of aerodynamic damping. The nominal aerodynamic damp- 

ing in torsion, Ex, is —0.0214, and the damping in plunging, E#, is 
0.8882 100 

5.3 Sensitivity of incidence angle to design variables, and resulting pertur- 
bations in stagger angle and reflex for a unit change in design variables. 101 

5.4 Computational times for present method using 129 x 49 node grid. . 106 

5.5 Change in steady flow quantities due to unit perturbations in ten design 
variables. The nominal steady lift, L, is 0.3942, the nominal drag, 

D, is —0.0168, the nominal moment about the leading edge, Mle, is 
—0.1837, and the nominal lift in the y-direction, Ly, is 0.3743 115 

5.6 Change in unsteady flow quantities due to unit perturbations in ten 

design variables. The nominal magnitude of the upstream pressure 
wave, |p up |, is 0.102, and the magnitude of the downstream pressure 
wave, |pdo W n|, is 0.249 118 

5.7 Computational times for present method using 129 x 49 node grid. . 125 


5.8 Change in steady and unsteady flow quantities due to unit perturba- 

tions in ten design variables near the hub of a modern fan exit guide 
vane. The nominal lift in the cascade direction, L y , is 0.4425, the nom- 
inal magnitude of the upstream pressure wave, |p up |, is 0.099, and the 
magnitude of the downstream pressure wave, |pdown|, is 0.350 129 

5.9 Change in steady and unsteady flow quantities due to unit perturba- 
tions in ten design variables near the midspan of a modern fan exit 
guide vane. The nominal lift in the cascade direction, L y , is 0.6579, 
the nominal magnitude of the upstream pressure wave, |p up |, is 0.209, 

and the magnitude of the downstream pressure wave, |pdown|, is 0.200. 131 

5.10 Change in steady and unsteady flow quantities due to unit perturba- 

tions in ten design variables near the tip of a modern fan exit guide 
vane. The nominal lift in the cascade direction, L y , is 0.8140, the nom- 
inal magnitude of the upstream pressure wave, |p up |, is 0.103, and the 
magnitude of the downstream pressure wave, |pdown|, is 0.263 132 


7 



Chapter 1 
Introduction 


1.1 Statement of the Problem 

As the efficiency of modern aircraft engines continues to increase, aeroacoustic and 
aeroelastic considerations will play an increasingly important role in the design of 
turbomachinery blading. Government regulations and community standards demand 
reduced levels of noise from aircraft, while competitive pressures require increased 
efficiency and reliability from modern designs. Unfortunately, aeroacoustic and aeroe- 
lastic performance can be very difficult to predict, due primarily to the complexity 
of the unsteady aerodynamic flowfield. In recent years, however, the computational 
modeling of these unsteady flows has substantially improved. 

Although current unsteady aerodynamic computational methods may improve the 
prediction of aeroelastic phenomena, they provide very little insight to the designer 
as to how to improve the aeroelastic behavior. As a result, the steady aerodynamic 
design and aeroelastic design phases during the development of fan, compressor, and 
turbine blading remain largely decoupled. After the airfoil shapes have been designed 
to satisfy their steady aerodynamic requirements, detailed aeroelastic studies are per- 
formed to determine whether the blades will meet design standards for flutter stability 
and fatigue. These studies can be very computationally expensive, particularly be- 
cause of the expense of determining the unsteady aerodynamic flowfield. If the blade 
fails to meet these requirements, the blade is redesigned, and the process is repeated. 
This redesign process increases the time and expense necessary to design a blade and 
misses an opportunity to design for steady and unsteady aerodynamic performance 
simultaneously. 

In current aeroacoustic analyses, the primary emphasis has been on the choice of 
the number of blades and vanes so that the so-called lower order modes do not prop- 
agate [1], Furthermore, steady blade loading, thickness, and camber effects are only 
approximated, if they are considered at all [2]. Hence, improved steady and unsteady 
aerodynamic modeling, particularly through the use of computational fluid dynamic 
(CFD) techniques, will lead to more accurate predictions of blade row response to 
incoming gusts, and help to illustrate the relationship between the blade shape and 
the radiated noise. 
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The goal of this research is to provide a framework for development and a real- 
istic implementation of useful design tools (as opposed to analysis tools) for design- 
ing turbomachinery blading for improved aeroacoustic and aeroelastic performance. 
These tools should be very computationally efficient, yet still model the dominant 
flow physics. In addition, these tools should provide physical insight that will lead to 
guidelines for future blade designs. 

1.1.1 Aeroelastic Problems in Turbomachines 

Aeroelastic phenomena arise from the interaction between inertial, elastic, and aero- 
dynamic forces [3]. For turbomachinery applications, this interaction may be illus- 
trated by examining a simplified version of the equation of motion for an airfoil 

TTIX -j- CX T kx — Fniotion(^') T Egu.stft') (1.1) 

On the left-hand side of Eq. (1-1), x is the displacement, the dots represent time 
derivatives, m is the mass of the blade, c represents the structural damping, and k is 
the blade stiffness. 

The right-hand side consists of two forces. The first force is due to self-induced 
oscillations and is a function of the displacement, velocity, and acceleration of the 
blade. If this force increases the energy of the vibration, the amplitude of the vibra- 
tion increases, and blade failure may occur. Such an unstable self-induced vibration 
is referred to as flutter. The second force on the right-hand side is due to an exter- 
nally excited oscillating motion (gust) where the force is independent of the blade 
displacement. The external forcing induces a response which in some cases may lead 
to high cycle fatigue failure of the blade, and is referred to as forced response [4]. 

These two aeroelastic phenomena, flutter and forced response, are two major 
types of aeroelastic problems encountered in modern aircraft engines. Both of these 
phenomena can lead to fatigue failure of one or more of the blades, and therefore are 
of great interest to aircraft engine designers. 

Flutter 

Self-excited blade vibrations that are sustained by the unsteady aerodynamic 
forces are of great concern to engine designers. If the aerodynamic forces produced 
by the blade vibration add energy to the blade motion, the vibration will grow expo- 
nentially until a limit cycle is reached or the blade fails. This aeroelastic phenomena 
is known as flutter. Flutter tends to occur near one of the lower natural frequencies 
of the blades. In modern engines, flutter may be encountered in a wide range of op- 
erating conditions. As a result, there are significant constraints on the design of fan, 
compressor, and turbine blading. Figure 1.1 shows the operating map for a typical fan 
or compressor, showing the boundaries for the most common types of flutter [3]. This 
figure illustrates the complexity of the design problem for a compressor. Deviations 
from the operating line can easily lead to one of the indicated types of flutter. Because 
the structural damping of the blades is usually small, the aerodynamic damping is of 
primary interest. The unsteady aerodynamics problem is to determine whether the 
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Figure 1.1: Axial compressor or fan characteristic map showing principal types of 
flutter and regions of occurrence (adapted from [3]). 

aerodynamic damping is positive (stabilizing) or negative (destabilizing) for a given 
blade motion and flow condition. If the system is unstable, then the blades must be 
redesigned. 

Forced Response 

Forced vibration is also a significant problem in the design of turbomachines. 
Forced vibrations may occur in fan, compressor, and turbine blades when a periodic 
aerodynamic forcing function acts on the blades in a given row. These forced vibra- 
tions may be destructive if the frequency of the forcing is near a resonant frequency of 
the blade. The forcing functions are generated at integer multiples of the engine rota- 
tional frequency and may arise from a number of different sources, the most common 
being aerodynamic interaction between adjacent blade rows. Specifically, unsteady 
aerodynamic forces often arise from the wakes of upstream blade rows impinging on 
the blade surface. Unlike the flutter problem, the aerodynamic forcing is independent 
of the blade motion, as indicated in Eq. (1.1). 

The frequencies at which these vibrations may occur can be predicted using a 
Campbell diagram [5]. Figure 1.2 shows a typical Campbell diagram for an axial-flow 
compressor, which displays the natural frequency of each blade vibration mode and 
the possible forcing function frequencies as functions of rotor speed. In this figure, 
the dashed lines indicate the natural frequencies of the first two bending and torsion 
vibratory modes of the rotor blades. Note that these curves increase somewhat with 
rotor speed due to centrifugal stiffening of the rotor. The solid lines are simply 
multiples of the rotor speed (so-called engine orders). It is at the crossing of these 
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Figure 1.2: Typical compressor resonance (or Campbell) diagram (adapted from [5]). 


two types of curves in the operating region that the possibility of destructive forced 
vibrations exists, because the forcing frequency corresponds to the resonant frequency 
of the blade. Hence, due to the small structural damping, large amplitude vibrations 
may be induced. 

Clearly, it is not possible to eliminate all sources of forced vibration from the 
operating range of a turbomachine, but designers need to either avoid these points 
or ensure that their elfect is minimal. For the forced response problem, an unsteady 
aerodynamic analysis is performed to determine the unsteady forces on the blades 
due to an incoming gust from the upstream blade row at the points where the curves 
cross on the Campbell diagram. If the unsteady aerodynamic forces induce large 
amplitude vibration at one or more of these points, the blade must be redesigned. 

1.1.2 Aeroacoustic Problems in Turbomachines 

Noise radiation continues to be an obstacle to the design and development of modern 
aircraft engines. Engine noise can be classified as either broadband noise (occurring 
over a continuous range of frequencies) or tonal noise (occurring at discrete frequen- 
cies). The relative magnitude of these two types of noise is primarily determined by 
the bypass ratio of the engine, i.e., the ratio of the amount of air flowing through 
the fan outside of the engine core to the amount of air flowing through the fan into 
the core. To increase the propulsive efficiency of engines, higher bypass ratios are 
desirable. As the bypass ratio of the engine increases, tonal noise becomes the more 
important contributor to the total engine noise signature [2], Furthermore, broad- 
band noise is a very complicated problem since it is usually caused by the interaction 
of random disturbances with blade surfaces and jet exhaust mixing. Because tonal 
noise is the larger contributor to the engine noise signature, broadband noise will not 
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be considered in this report. 

The source of tonal noise is nonuniformities in the flowfield such as inlet dis- 
tortion, inlet angle of attack, or convected wakes from upstream blade rows. The 
noise is produced by an interaction between the flow nonuniformity (or “gust”), the 
blade row, and the duct surrounding the blade row coupling to produce propagating 
acoustic modes. Currently, most tonal noise reduction methods are aimed at either 
reducing the strength of the gust by increasing the distance between blade rows, or 
by modifying the duct through the use of casing treatment. In addition, blade counts 
are carefully chosen. Pressure waves that decay as they propagate are referred to 
as “cut-oif.” Propagating pressure waves are referred to as “cut-on,” and it is these 
waves which result in audible noise. Acousticians use analytical models of the duct to 
choose blade and vane counts to prevent lower order pressure modes from propagat- 
ing [1]. All three of these noise reduction techniques, increasing blade row spacing, 
casing treatment, and modifying blade/ vane counts are not desirable because these 
methods add weight and expense to the finished engine. 

The unsteady aerodynamics problem is to determine the response of a blade row 
to an incoming gust, determine the. pressure field generated by the wake-blade in- 
teraction and determine whether these pressure waves will propagate outside of the 
engine. Most unsteady aerodynamic methods, however, are best viewed as analysis 
tools rather than design tools. They are capable of solving the direct problem where 
the shape of the airfoil as well as the flow conditions are specified. Unfortunately, 
except through trial and error or extensive parametric studies, these analyses do not 
provide physical insight into how, for example, to design blade shapes to reduce their 
acoustic response without compromising their aerodynamic efficiency. An additional 
complication is that the frequencies associated with aeroacoustic analysis are ap- 
proximately an order of magnitude higher than those found in aeroelastic problems, 
requiring finer discretization of the governing equations of the unsteady flowfield and 
increased computational expense. 


1.2 Previous Work 

1.2.1 Unsteady Aerodynamic Models 

The first unsteady aerodynamic analyses were semi-analytical methods that were 
designed to determine the unsteady surface pressure on a cascade of two-dimensional 
flat plates due to bending or torsional vibration. In the semi-analytical approach, 
a number of point singularities are placed on the blade surface, and an appropriate 
kernel function is used to calculate the unsteady flowfield. One of the first such 
analyses was performed by Whitehead [6], who considered uniform, incompressible, 
undeflected flow over the plates. Because the flow is not deflected, however, there is no 
steady loading on the blade, and the model failed to predict experimentally observed 
bending flutter. To include the effect of steady loading, Whitehead later allowed the 
blades to deflect the mean flow, resulting in steady loading on the blades [7], and 
bending flutter was predicted. This provided numerical evidence that steady blade 
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loading (and therefore the blade shape design) had an important effect on flutter 
in aircraft engines. Later, Atassi and Akai [8] distributed point singularities along 
the surface of an airfoil with finite camber and thickness. This model allowed them 
to analyze incompressible flows through more realistic cascades, also showing the 
importance of steady loading and blade shape on flutter. 

Both Whitehead [9] and Smith [10] investigated subsonic compressible flows over 
flat plate airfoils. Smith’s approach was developed primarily for noise analysis in an 
attempt to determine the propagation characteristics of acoustic modes generated by 
a cascade of airfoils. Smith also measured the acoustic waves experimentally, and 
found that the theory and experiment agreed well for the unloaded blade case, but 
the theory was inadequate to deal with the extra sound sources introduced by steady 
blade loading, and underestimated the generated wave amplitude in this case. Both 
Whitehead and Smith found that steady loading is extremely important to predict 
correctly the amplitudes of the outgoing pressure waves. 

Advanced designs required researchers to analyze transonic and supersonic flows. 
Adamczyk and Goldstein [11] and others [12, 13, 14] investigated vibrating flat plates 
in supersonic flow which is axially subsonic. Unfortunately, these models did not 
predict flutter at the experimentally observed frequencies and Mach numbers. More 
sophisticated models, such as the in-passage shock model developed by Goldstein, 
Braun, and Adamczyk [15], analyzed supersonic cascades with finite strength shocks. 
The in-passage shock model predicted experimentally observed supersonic bending 
flutter, at least for low reduced frequencies. The design problem was also beginning 
to be considered. Bendiksen [16] used a perturbation analysis to include effects due 
to steady loading, thickness, camber, angle of attack, and shock motion, showing that 
these effects are important in flutter prediction. 

Finally, some three-dimensional problems have been considered using a semi- 
analytical approach. Namba [17] has developed a lifting surface analysis to determine 
the unsteady flow over vibrating three-dimensional flat plate cascades in the sub- 
sonic, transonic, and supersonic regimes. Although his approach demonstrates the 
importance of three-dimensional effects, its application is primarily limited to lightly 
loaded fan blades. 

Although these semi-analytical methods have developed greatly over the years, 
they are insufficient to analyze most unsteady aerodynamic problems in modern en- 
gines, especially at off-design conditions. The blades in modern designs are heavily 
loaded, and violate a number of the assumptions used in the semi-analytical approach. 

As computers became more powerful, so-called field methods became an impor- 
tant research area. Using this approach, a set of governing equations (i.e., potential, 
Euler, or Navier- Stokes) is solved using computational fluid dynamics techniques. 
This approach has the advantage that many of the effects not included in the analyt- 
ical models may be easily incorporated, e.g., arbitrary blade geometries, complicated 
shock structures, and various flow models. 

In the area of field methods, two main approaches have emerged. These are re- 
ferred to as time-marching and time-linearized methods. Using the time-marching 
approach, the unsteady governing equations of the fluid motion are time-accurately 
marched subject to some appropriate set of unsteady boundary conditions. For ex- 
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ample, to analyze a flutter problem, the airfoils are prescribed to vibrate at some 
fixed frequency and interblade phase angle. Once the initial transients have decayed, 
the unsteady flowfield is found to be periodic in time. The advantage to this ap- 
proach is that nonlinear flow effects are incorporated, and finite blade motions may 
be considered. 

Two-dimensional, time- accurate Euler solution methods have been developed by 
Giles [18] and others. Three-dimensional unsteady flow problems have been in- 
vestigated by Ni and Sharma [19] using a time-accurate Euler method, but their 
results require several hours of supercomputer time to obtain a solution. Other 
three-dimensional methods have been developed to investigate rotor-stator interac- 
tion. Saxer and Giles [20] used an Euler method and Rai [21] used a Navier-Stokes 
analysis. Although these methods model many of the important physical mechanisms 
of the unsteady flow problem, the massive amount of time required to march these 
equations time-accurately will prevent this type of analysis from being used for design 
purposes for many years to come, especially for three-dimensional problems. 

The second approach in field methods is the use of linearized analyses to investigate 
small perturbation unsteady flows. Using this approach, the flow is assumed to be 
composed of a nonlinear mean or steady flow plus a small unsteady perturbation flow. 
The linearized equations which describe the unsteady perturbation are linear variable 
coefficient equations in the unknown complex amplitude of the harmonic motion of 
the flow. 

Verdon and Caspar [22] determined the mean flow through a subsonic cascade 
using a steady full potential method. They then linearized the potential equation 
about the mean flow to solve for the small disturbance unsteady flow. Whitehead 
and Grant [23] and Hall [24] used a similar approach but used finite elements instead 
of finite differences to discretize the linearized equations. Verdon and Caspar later 
extended their model to transonic flows using shock fitting to model the shock mo- 
tion [25]. Although these potential methods work well for subsonic flutter problems, 
their use for forced response and aeroacoustic applications has been limited because 
the model does not in general include unsteady vorticity, as would have to be included 
to determine the response of a blade row to an incoming vortical gust. 

Within a linearized potential framework, Goldstein [26] developed a velocity split- 
ting technique that includes unsteady vorticity and entropy. Goldstein assumed an 
isentropic and irrotational mean flow, and split the unsteady velocity into rotational 
and irrotational parts. The rotational unsteady velocity, which contains the unsteady 
vorticity, is expressed analytically using the drift and stream functions from the steady 
flow solution. The rotational velocity appears in an inhomogeneous source term in 
the linearized potential equation. Hence, unsteady flows with vorticity could be an- 
alyzed without having to numerically solve the Euler equations. Unfortunately, this 
analysis could only be applied to flat plate airfoils. Atassi and Grzedzinski [27] de- 
veloped a modification to this method that allowed real airfoils to be analyzed. Hall 
and Verdon [28] implemented this technique and showed that vortical gust effects 
could be modeled accurately without the computational expense of solving the Euler 
equations. 

Even with vortical gust extensions, the linearized potential formulation has its 



1.2. PREVIOUS WORK 


15 


limits. For transonic flows with strong shocks, the Euler equations must be solved. 
Also, the potential formulation does not translate well to three-dimensional applica- 
tions due to the assumption of irrotational steady flow. For these and other reasons, 
linearized methods continue to advance. 

Hall and Crawley [29] originally solved the linearized Euler equations using a 
Newton iteration technique to solve the nonlinear steady Euler equations and a direct 
method to solve the unsteady linearized Euler equations. This approach also used 
shock fitting to model the shock motion. There were some problems with this method, 
however. For example, the direct solution method prevented fine computational grids 
from being used on typical workstation computers, particularly if the analysis were 
to be extended to three dimensions. To circumvent this and other problems, Hall 
and Clark [30] and Holmes and Chuang [31] implemented linearized Euler methods 
using an iterative solution method rather than a direct solution. In addition, Hall, 
Clark, and Lorence [32] and Lindquist and Giles [33] showed that transonic aeroelastic 
problems can be modeled appropriately within a linearized framework using shock 
capturing. Using a similar technique, Clark and Hall [34] have implemented a method 
for solving the linearized Navier-Stokes equations so that stall flutter problems may 
be analyzed. Finally, Hall and Lorence [35] extended the linearized Euler analysis to 
three dimensions and showed the importance of three-dimensionality on the unsteady 
flow in fans. 

Linearized methods require significantly less computational time than their time- 
marching counterparts. For most unsteady aerodynamic flows of interest, the lin- 
earized approach requires one to two orders of magnitude less computational time 
than an equivalent time-marching calculation. Although for some complex appli- 
cations, time-accurate time-marching methods may be employed, in most cases, a 
linearized analysis is sufficient to model the unsteady flowfield. 

1.2.2 Design Methods and Optimization 

A substantial body of work exists on the inverse design and optimal design of airfoils. 
Most of this work, however, is directed at achieving desirable steady flow properties. 
In an inverse method, a pressure distribution is specified and the analysis produces 
an airfoil shape (if possible) that will produce the desired distribution. For exam- 
ple, Lighthill [36] developed an inverse design method based on conformal mapping 
techniques. More recently, a number of investigators have proposed inverse design 
techniques based on modern computational fluid dynamic algorithms (e.g., [37]). 

More popular are optimal design techniques, where an initial airfoil shape is spec- 
ified, and the analysis attempts to minimize some quantity (e.g., steady aerodynamic 
losses) while satisfying some appropriate constraints (e.g., the blade maintains the 
correct turning). A number of investigators have used nonlinear programming tech- 
niques to solve this problem (e.g., [38]), and Jameson has suggested that this problem 
may be viewed as an optimal control problem [39]. 

Researchers have also developed some optimization techniques for aeroelastic 
problems in turbomachinery. For example, Crawley and Hall [40] developed a method 
to calculate an optimal distribution of “mistuning” for a blade row. These techniques, 
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however, have focused on structural optimization rather than optimization of the un- 
steady aerodynamic behavior. 

One of the key ingredients in optimization algorithms is the evaluation of the 
sensitivity of the quantity to be optimized (for example, the flutter stability or ef- 
ficiency of a cascade) to a small change in a physical parameter (such as the airfoil 
shape). Sensitivity analysis of structures has been an active area of research for the 
past decade [41, 42], Recently, researchers have begun to develop similar sensitivity 
analysis techniques for steady aerodynamic problems. For example, Taylor et al [43] 
and Baysal and Eleshaky [44] have computed the effect of modifying the shape of a 
nozzle on the flow in the nozzle. Their work was based on a sensitivity analysis of 
the discretized Euler equations. Most recently, such techniques have been applied to 
airfoil design [45]. Despite these advances, only a few unsteady sensitivity analyses 
have been reported, for example the semi-analytical panel method of Murthy and 
Kaza [46], Other unsteady sensitivity analyses have been computed by differencing 
two slightly different nominal unsteady solutions. The use of these types of finite 
differences, as opposed to analytical perturbations, are not as desirable because of 
their increased computational expense and susceptibility to truncation error. 


1.3 Present Method 

In this report, strategies for efficient aeroacoustic and aeroelastic design of turboma- 
chinery blades will be examined. An essential part of these strategies is a new method 
for computing the sensitivity of unsteady flows in cascades due to small changes in 
blade geometry. The approach is general in nature and may be applied to different 
governing equations and numerical schemes. Consequently, an appropriate equation 
set must be chosen. For many aeroacoustic and aeroelastic problems in turbomachin- 
ery, an inviscid flow analysis is sufficient to model the dominant physical mechanisms 
of the unsteady flowfield. Hence, a linearized potential or linearized Euler framework 
would be applicable. Since vortical gusts must be analyzed, however, the traditional 
linearized potential method is not appropriate. The linearized Euler approach, while 
incorporating complete inviscid physics, is still computationally expensive to solve 
directly, as the present method requires. Therefore, a linearized potential model with 
vortical gust extensions is the logical choice. This is the method that will be used in 
this report to demonstrate the sensitivity analysis procedure. 

Specifically, the nominal analysis used here is a linearized harmonic unsteady po- 
tential method based on a deforming grid variational principle and finite element 
method by Hall [24] extended using rapid distortion theory to include the effect of 
incident vortical gusts due to wake interaction. The sensitivity analysis procedure is 
as follows. First, the nominal flow equations are discretized and solved using a finite 
element procedure. Next, a perturbation analysis is performed on the discretized 
equations from the nominal finite element analysis. This leads to a set of linear equa- 
tions for the sensitivity of the unsteady potential due to small changes in the airfoil 
shape. If the nominal unsteady analysis is computed with LU decomposition, the sen- 
sitivities may then be computed by back-substitution. Consequently, the sensitivity 
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of the unsteady potential may be computed very efficiently. 

Once the sensitivities have been calculated, a range of approaches exist for improv- 
ing the aeroacoustic and aeroelastic performance of the cascade. For example, the 
calculated sensitivities may be used by themselves as a guide to redesign the airfoil. 
A more thorough approach would be to perform an optimization on the aeroacous- 
tic or aeroelastic behavior of the cascade, subject to appropriate constraints. The 
radiated noise, for example, could be minimized subject to the constraint that the 
desired steady flow turning is achieved. The results shown in this report will essen- 
tially be a compromise between these two methods, where the calculated sensitivities 
are post-processed such that design changes satisfy steady aerodynamic constraints. 
The framework for the use of this technique in more complicated design situations 
will also be discussed. 


1.4 Outline of Report 

fn Chapter 2, the governing equations and basic theory of the nominal flow analysis 
will be presented. This includes both the steady flow potential equation and the 
linearized potential equation and the vortical gust extension. Use of a deforming 
computational grid will be discussed, as will the theoretical description of the steady 
and unsteady boundary conditions. 

Chapter 3 contains a discussion of the basic numerical discretization scheme, based 
on Hall’s finite element method. The grid generation procedure will be discussed, as 
will the implementation of the boundary conditions, with particular emphasis on the 
nonreflecting far-field boundary conditions. Finally, the matrix assembly and solution 
will be described. 

In Chapter 4, the sensitivity analysis procedure will be examined in detail for both 
the steady and unsteady flow equations. Although the procedure will be applied to 
the present potential method, this chapter should clarify how the sensitivity analysis 
would be performed on other governing equations and discretization schemes. 

In Chapter 5, results of the sensitivity analysis will be presented, demonstrating 
how the analysis may be used to redesign blade rows for improved aeroelastic and 
aeroacoustic performance. Also, the computational efficiency of the method will be 
discussed, as will its range of effectiveness. 

In Chapter 6, three-dimensional flow considerations will be examined. Specifical- 
ly, application of the analysis in an actual design environment will be discussed , with 
particular emphasis on how three-dimensional problems could be analyzed using the 
present method. 

Finally, in Chapter 7, some conclusions from the present analysis will be presented 
as well as some recommendations for future work in this area. 



Chapter 2 

Nominal Flow Field Description 


In this chapter, the equations and boundary conditions governing the steady and 
small disturbance unsteady flow through a two-dimensional cascade of airfoils are 
introduced. Section 2.1 contains a description of the rapid distortion theory used 
so that unsteady vortical flows may be modeled within a potential flow framework. 
This theory results in a set of sequentially coupled partial differential equations that 
must be solved numerically to obtain the unsteady flow. In Section 2.2, an extension 
of a variational principle originally developed by Bateman [47] is described. This 
variational principle has as its Euler-Lagrange equation the potential flow equation 
developed in Section 2.1. This variational principle will be used in Chapter 3 to 
construct a finite element description of the steady and unsteady flow. Finally, in 
Sections 2.3 and 2.4, the steady flow and small disturbance unsteady flow boundary 
conditions will be discussed, including the far-field nonreflecting boundary conditions. 


2.1 Euler Equations 


The equations governing the fluid motion through a cascade of airfoils may be derived 
from the integral conservation laws of mass, momentum, and energy, along with the 
equation of state for an ideal gas. For this analysis, the flow is assumed to be inviscid 
and compressible, and the fluid is assumed to be a perfect gas with specific heat ratio 
7 = Cp/cy (a complete list of symbols used in this report is given in Appendix A). 
Consequently, the governing equations of the fluid are the Euler equations, which in 
conservation form are 

g + V-(pv) = 0 


dt 


4- V • pvv -f Vp = 0 


( 2 . 1 ) 
( 2 . 2 ) 
(2.3) 

where p is the density, v is the velocity, and p is the pressure. These quantities are 
related by the state equation for a perfect gas 



p = pRuTf 


(2.4) 
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where EQ is the universal gas constant and Tf is the temperature. The superscript 
indicates that the quantity is unsteady. For this discussion, it is useful to define the 
entropy of the fluid, sj. This is most easily accomplished using the thermodynamic 
relation defined by Gibbs, 

Tf dsj = Cp dfj - \ dp (2.5) 

P 

The temperature dependence may be removed through the use of the state equation, 
Eq. (2.4), resulting in 

d P d P 

dsf-c v —-c p — (2.6) 

VP 

Equation (2.6) may then be integrated to obtain the entropy change between two 
fluid states. 

It is further assumed that the nonlinear time-varying flow may be split into two 
parts: a nonlinear mean or steady flow and a small unsteady perturbation flow that 
is harmonic in time. In other words, each of the variables in the above equations may 
be expanded in a perturbation series. Furthermore, it is assumed that the steady flow 
is irrotational and homentropic, so that the steady velocity may be written as the 
gradient of a scalar velocity potential, <fr. Under these assumptions, the perturbation 


series for the flow variables may 

be written as 



p{*>y,t) = 

R(x,y ) 

+ 

p(x,y)e JU,t 

v(x,y,t) = 

V$(x,y) 

+ 

v(x,t/)e JU,£ 

II 

>5 

H 

P{x,y) 

+ 

p(x,y)e Ju;t 

s f {x,y,t) = 



s f (x,y)e 3u,t 

' 

total flow 

nonlinear 


small harmonic 


(2.7) 


mean 


disturbance 


To obtain the steady and unsteady flow equations, these expansions are substi- 
tuted into the governing equations, Eqs. (2. 1)— (2.6) , collecting terms of equal order. 
The zeroth-order terms result in the steady flow equations, while the first-order terms 
describe the unsteady flow. The zeroth-order continuity equation is 


V-(f?V$) = 0 


(2,8) 


The zeroth-order momentum equation and equation of state may be integrated 
and combined to obtain the steady form of Bernoulli’s Equation, i.e., 


R = pt 




i-i 


(2.9) 


or 


P =PT 


1 - 


i-i 


(2,10) 
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where pj and pr are the total pressure and density, respectively, and Ct is the total 
speed of sound. Substituting Eq. (2.9) into Eq. (2.8) results in an equation for the 
mean flow that is only a function of the potential and the steady speed of sound, 


V 2 $ = — 


C 2 L2 


-V$ • V(V$) 2 


( 2 . 11 ) 


where 

C 2 = C\ - (V$) 2 (2.12) 

Note that the steady flow is completely described by a single equation for the scalar 
potential $. It is indeed nonlinear, as was noted earlier. 

Having described the equation governing the steady flow, we now consider the 
unsteady small disturbance equations. Collection of the first-order terms results in the 
following set of equations describing the behavior of the small disturbance unsteady 
flow 


(v - s } V<I>/2) + [(v - syV$/2) •V]VHVQ=0 (2.14) 

Dt {rCp) + j R 7 ' ^ = ° ( 2 ' 15 ) 

where D/Dt = d/dt + V$ • V is the convective derivative. Note that while the steady 
flow may be described using a single equation, the unsteady flow is governed by a 
system of four simultaneously coupled equations. Computationally, this description 
of the unsteady flow would result in a very large system of equations to solve. Hence, 
a more compact description of the unsteady flow is desired. 


Dsf 

rw 


= 0 


(2.13) 


2.1.1 Goldstein Decomposition 


Equations (2.14) and (2.15) may be simplified further using the Goldstein velocity 
decomposition [26]. The unsteady velocity, v, is split into a rotational part, v R , and 
an irrotational part that is written as the gradient of a scalar velocity potential, <f>, 
so that 


v = + v H 


(2.16) 


It should be noted that there is no unique choice for v R and </>, since it is the combi- 
nation of the two that results in the actual unsteady flow velocity. However, if v R is 
chosen such that 


D R 
—v + 
Dt 


,R 


v ] v$ + v (i ) =0 


(2.17) 


then the unsteady pressure is only a function of the unsteady potential, so that 
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Substituting Eqs. (2.16) and (2.18) into the linearized momentum and continuity 
equations [Eqs. (2.14) and (2.15)] results in a system of equations for the unsteady 
entropy, rotational velocity, and velocity potential, i.e. , 


Ds j 
~Dt 


= 0 


(2.19) 


A (y R -s f Vt>/ 2) + [(v R — 5/V<&/2) -V] V# = 0 (2.20) 

Note that the equations for entropy, rotational velocity, and velocity potential are now 
only sequentially coupled. Although four equations still must be solved, they need 
not be solved simultaneously. As a result, a numerical solution to these equations 
may be obtained for considerably less computational cost than the earlier set of 
simultaneously coupled equations. Furthermore, in this report, the effect of unsteady 
entropy will not be investigated, so it is assumed that sj = 0. Consequently, only 
Eqs. (2.20) and (2.21) must be solved. 

Goldstein [26] showed that Eq. (2.20) may be solved analytically for certain cases 
in terms of two functions. One is the well-known stream function, \l>(x, y). The other 
is the drift function, A (x,y), which is defined as the time it takes a mean flow fluid 
particle to move between two points on a streamline, i.e., 


A (x,y) = A(x 0 , yo) + 


L 


(*>v) 1 


(®o,vo) |V$| 


difr 


( 2 . 22 ) 


where the location (xo, yo) is a point upstream of (x, y) on the same steady streamline, 
|V$| is the magnitude of the steady flow velocity, and the dilferential distance dijj is 
measured along the streamline. The drift and stream functions are often referred to 
as the Lagrangian coordinates of the fluid. 

Figure 2.1 shows contours of the drift and stream functions for a typical fan 
exit guide vane. Note that upstream of the cascade, the drift function contours are 
aligned between the blade passages, while downstream of the airfoils the contours are 
no longer aligned. This indicates that there is circulation around the airfoil, i.e., the 
steady velocities are different on the two sides of the airfoil surface. Furthermore, 
note that the drift contours are well-behaved in the interior of the blade passages, 
but change rapidly near the blade and wake surfaces. This is due to the presence of a 
stagnation point at the leading edge of the airfoil. Examining Eq. (2.22) shows that if 
the steady velocity is zero at some point (i.e., a stagnation point), the drift function 
becomes infinite. This property of the drift function is problematic in the derivation 
of a general expression for the rotational velocity over airfoils, so a detailed analysis 
of the drift function behavior near a stagnation point is warranted. 

Consider the steady flow field near a stagnation point. Sufficiently close to a solid 
boundary, the boundary appears to be a plane surface (unless the boundary happens 
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Figure 2.1: Contours of drift function, A(x,y), and stream function, ty(x ,y), for a 
typical fan exit guide vane. 


to have a discontinuity of slope there) [48]. If the flow is incompressible, the governing 
equation for the steady potential is Laplace’s equation, i.e., 


V 2 $ = 0 (2.23) 

The solution that satisfies the boundary conditions is 

4 = \m,x 2 - y 2 ) (2.24) 

where A is some constant. Alternatively, the flow may be described in terms of the 
stream function by 

^ = Axy (2.25) 

Figure 2.2 shows contours of the drift and stream functions near a stagnation point. 
The stream function is shown near a plane boundary [i.e., a flow described by Eq. (2.25)], 
while the drift function is shown near the stagnation point of an actual airfoil. 
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Figure 2.2: Left, contours of the stream function, near a solid plane boundary. 
Right, contours of the drift function, A, near an airfoil stagnation point. 


Returning to the flow near a plane boundary, it is clear from Eq. (2.24) that the 
velocity in the direction normal to the surface (the x-direction) is 

d§ 

S^ = Ax < 2 ' 26 ) 

So the velocity is proportional to the distance from the boundary. Hence, using the 
definition of the drift function [Eq. (2.22)], near the surface 


A = — - In x 

A 


(2.27) 


So there is a logarithmic singularity in the drift function near the stagnation point. 
In general, then, we may define a constant geo such that 

/aiviiy 1 


Uq — ~ ’ 


{ dn ) 


(2.28) 


sp 


where |V$| is the magnitude of the steady velocity, n is measured normal to the 
airfoil surface, and the derivative is evaluated at the stagnation point. Using this 
definition, from Eq. (2.27) it is apparent that to leading order 


3A _ ao 
dn n 


(2.29) 


Furthermore, using Eq. (2.25), Eq. (2.29) may be rewritten in terms of the stream 
function, so that 

<9A o 0 


d$ # - # 0 


+ ... 


(2.30) 
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where is the value of the stream function on the stagnation streamline. This result 
will prove to be useful shortly when a general expression for v R will be derived. First, 
however, we will return to Goldstein’s description of the rotational velocity. 

Goldstein’s paper did not explicitly address the singular nature of the drift func- 
tion, so we will put this property aside for the moment. In general, the Goldstein 
rotational velocity, v G , may be written in terms of the drift and stream functions as 

v G = ( Cl VA + c 2 V'5)exp[j'(/GA + I<^)\ (2.31) 

Here K 1 and K 2 are essentially wave numbers along drift and stream contours, and 
Ci and c 2 are constants. Goldstein chose Ci and c 2 so that the rotational velocity in 
the upstream far field is divergence-free, i.e., 

V ■ v G = 0 (2.32) 

or 

c\K\ + c 2 K 2 = 0 (2.33) 

The constants c\ and c 2 may then be uniquely determined by specifying the magnitude 
and phase of the vorticity of the incoming gust. 

With v G defined, in principle all that remains is to solve Eq. (2.21) for the unsteady 
potential, <j). However, for flows over bodies that contain a stagnation point, this 
formulation leads to a singular rotational velocity on the blade and wake surfaces due 
to the logarithmic singularity in the drift function at the stagnation point described 
earlier. To avoid this singularity, one possible approach would be to set c a to zero in 
Eq. (2.31) to remove the terms that depend on the gradient of the drift function, and 
choose c 2 to match the specified vorticity. Unfortunately, although this modification 
produces a bounded rotational velocity, the rotational velocity has indeterminate 
phase along the blade and wake surfaces. As a result, the divergence of the rotational 
velocity still has a strong singularity at the airfoil and wake surfaces. Since the 
source term in Eq. (2.21) is dependent on the divergence of the rotational velocity, a 
modification to the description of the rotational velocity is necessary. 

2.1.2 Atassi Decomposition 

Atassi and Grzedzinski [27] developed a modified velocity splitting that is uniformly 
valid for flows around bodies. Shortly thereafter, Hall and Verdon [28] used this 
technique to study unsteady flows in cascades. 

Atassi and Grzedzinski suggested adding the gradient of a convected, and therefore 
pressureless, scalar potential, to Goldstein’s rotational velocity, so that 

v R = v G + V4> (2.34) 

Since the curl of the gradient of a scalar is zero, this additional term does not change 
the unsteady vorticity in the flow. 
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This potential may not be simply chosen arbitrarily. First of all, this new form of 
v R must satisfy Eq. (2.20). This is true if 



(2.35) 


which implies that it has no pressure associated with it [see Eq. (2.18)]. Since it is 
convected, the potential propagates the same way as the Goldstein rotational velocity, 

i- e -, 

<£ = $(^r)exp[j(/GA + /F 2 ^)] (2.36) 

where $ is an as yet undetermined function. Using the form of <j> given in Eq. (2.36), 
the rotational velocity may be written as 


v R = 


Cl + 


dct> 

dA 


VA + 




ex p[;(A'i A + /{»!■)] 


(2.37) 


We wish to determine the potential (f> such that v R is zero on the blade and 
wake surfaces, which will produce the desired result that the normal and tangential 
components of the rotational velocity field are regular along the blade and wake 
surfaces. In terms of the unit surface tangent, s, and the unit surface normal, n, we 
want 


v fi • s = (v G + • s = 0 


(2.38) 


and 

v R • n = (v G + V<^) • n = 0 


(2.39) 


The first condition, Eq. (2.38), is satisfied if $ is a constant. To illustrate this, 
consider the tangential component of v R . Using Eq. (2.37), if $ is constant, then 


/ . ~\ dA I ~\ chk 

( Cl +]K,$) - 7 - + (c 2 + jK 2 <$>) — = 0 


(2.40) 


Since the airfoil surface and wake is a stagnation streamline, d'h/ds is zero. Hence, 
if $ is chosen such that 

* = f; ( 2 - 41 ) 

then Eq. (2.40) [and therefore Eq. (2.38)] is satisfied. This choice of $ removes the 
dominant singularity in v G which is on the order of VA. In fact, this is equivalent 
to setting ci to zero in Eq. (2.31) and choosing c 2 to match the specified vorticity, as 
suggested earlier. 

Unfortunately, v R may still have indeterminate phase. The singular behavior of 
the normal component of v R is only completely eliminated when <j> satisfies Eq. (2.39). 
To eliminate this singular behavior, we will proceed by splitting <^> into two parts, so 
that 

4> — 4>\ A 4>2 (2.42) 

where, fa is the potential associated with the constant <1 given in Eq. (2.41), and <f> 2 
is still to be determined. So as not to violate Eq. (2.38), ^ 2 will be chosen so that its 
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value will be zero on the blade and wake surfaces, while allowing <j 2 to be dependent 
on the stream function, '5. 

Considering both parts of fa the normal component of v R may be written as 



dfa dfa \ 
dA + dA J 




d fa 
d'F + 




(2.43) 


In light of the choice 


of fa, this reduces to 


dfadA _ ( c,/{ 2 , 3$ n 

dA dn + \ C2 I< x + d'H ) dn ~ ° 


(2.44) 


or 


dfa dA , c x I< 2 , dfa\ d'F n 

dA d < t + C2 A, + a® ) dn ~ ° 


(2.45) 


Next, we wish to 
surface of the airfoil, 
point gives 


determine the behavior of the first term in Eq. (2.45) near the 
Performing a Taylor expansion of dfa/d A about the stagnation 


dfa 

dA 


d 

dV 



(tf - tf„) + ... 


(2.46) 


This equation, combined with the expression for dAjdlt obtained in the earlier dis- 
cussion of the behavior of the drift function near a stagnation point [Eq. (2.30)], leads 
to 


' d 

( dfa) 

, c x I< 2 

dfa 

dV 

dF 

\dA) 

a o + c 2 T , 

d^ 

dn 


(2.47) 


As was shown earlier, fa has the form 


fa ~ $ 2 (^)exp[j(A' 1 A + K 2 ty)\ 
Substituting this expression into Eq. (2.47) gives 


(1 + ja 0 K x )^ + C2 -^+ jK 2 * 2 


d$ 


K ! 


dH A 

dn ~° 


(2.48) 


Also, because the cascade is periodic, a similar boundary condition must be applied 
on the next adjacent blade surface. Hence, the value of the stream function at this 
surface is required. The stream function may be computed by integrating the mass 
flux over the blade-to-blade gap 


\?(C?) — tfo + 



(2.49) 


where G is the blade-to-blade gap measured in the y-direction, and Q is the mass 
flux. Typically, the upstream steady density, f?_oo, the upstream free stream velocity, 
K-oo, and the upstream flow angle, <*,, are specified. Using these values, the stream 
function may be written as 


#( X,y + G) = #0 + R-ooV-ooG cosf2_oo 


(2.50) 
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Consequently, since the potential should vanish at the blade and wake surfaces, we 
let $ 2 take the form 


$ 2 ('^) = B sin 


2tt(^ - $p) 
R—oo^^—coCr COS fl — oo 


(2.51) 


where B is a constant independent of \I/. Substituting this functional form into the 
boundary condition [Eq. (2.48)] gives 


2xB(l + japKi) 

R. — QQ f qq G COS fl — r 


COS 


2tt($ _ $„) 


+ c 2 - + jK 2 B sin 


solving for B gives 


K\ 


B = 


R. QQ k-QQ Cj COS H | 

27 r(tf - %) 


J 'Ps'Po 


R— OO k— oqC COS Oj — r 


= 0 


'P='I , 0 


R—ooV—qqG COS 0 — oo(c 2 R\ Ci K 2 ') 


2-kI<i{ 1 + ja 0 Ki) 

After performing some algebra, the final expression for ^ is obtained, i.e. , 




jR — CO V-ooGcosCl-oofaKi - c x I< 2 ) . 

~ , 7 sin 

27 r(l + ja 0 I<i) 


2tt(^ - ^q) 


(2.52) 


(2.53) 


R — (V) V-co G COS fl — r 


x exp[j(Ai A + K 2 ty)] (2.54) 

This expression for <f> results in a rotational velocity that is regular along the blade 
and wake surfaces, so that the divergence of the rotational velocity may be calculated 
despite the presence of a stagnation point. Note that the complete expression for v R 
is relatively simple, and may be calculated once the steady flow has been computed 
and the desired magnitude and phase of the unsteady vorticity has been chosen. Now 
that a uniformly valid expression for v R has been written, then, the next task is to 
calculate the unsteady velocity potential, using Eq. (2.21). 


2.2 Extension of Bateman’s Variational Principle 

Now that the governing equations have been developed, it is useful to consider how 
these equations will be solved numerically. Essentially there are two choices: finite 
differences and finite elements. Both of these methods have been used to solve the 
linearized potential equation for blade motion and pressure gust analyses [22, 23]. 
In addition, Hall and Verdon used finite differences to solve the linearized potential 
equation with rapid distortion theory to analyze vortical gusts [28]. For this report, 
a finite element discretization of the field equations is used. There are three main 
reasons for this approach. First, finite element techniques are versatile, elegant, and 
relatively easy to implement. Second, it is believed that a finite element approach will 
reduce the overall truncation error because the governing equation is solved in the 
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integral sense instead of the differential sense. Truncation errors tend to be higher 
in derivative evaluation than integral evaluation. Finally, a potential method for 
blade motion analyses developed by Hall [24] was readily available. The details of 
the finite element procedure will be given in Chapter 3. This section will describe the 
variational principle on which the finite element method is based and its application 
to the current work, including the effect of unsteady vorticity. 


2.2.1 Nonlinear Full Variational Principle 

First, consider irrotational and homentropic flow. Bateman [47] developed a varia- 
tional principle that states that for a steady flow, the volume integral of the pressure 
must have an extreme value, and the associated Euler- Lagrange equation is the steady 
conservation of mass. Hall [24] later extended this principle for application to cas- 
cades by considering temporally periodic flow. In this case, the pressure integrated 
over a domain E and over a period A is extremized. In functional form, 


n 


=i/ rt 

A Ja J is 


p dx dy dt + 




ds dt 


(2.55) 


where Q is the prescribed mass flux on the boundary and s is the distance along the 
boundary T. Taking the first variation of n and setting to zero gives 


<!>n = J J J 8p dx dy dt + J <j> Q 8$ ds dt ~ 0 

Bernoulli’s Equation says that 


(2.56) 


+ Tt- C 


(2.57) 


where C is a constant. Since we have assumed that the flow is homentropic, the 
variation of the pressure may be written as 


/ . * d 0 

8p = —p ( • S78(f) -)- — <*><^ 


(2.58) 


Substituting this expression into the equation for (ffl, using the divergence theorem 
and integration by parts gives 

{n = 1 i SL [ v ' WI) + 1 ] s * di d y dt 


+ 6idsdt = 0 


(2.59) 


To extremize n, dn must be zero for all permissible variations in q l. Therefore, 
Eq. (2.59) says that the conservation of mass must be satisfied in the domain E. On 
the boundary T, Dirichlet boundary conditions may be imposed, since then 8<j> will be 
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zero, or the Neumann boundary condition may be used, so that pd^jdn = Q. This 
is the so-called natural boundary condition. 

Since the steady and unsteady flow is governed by the conservation of mass, this 
variational principle is the appropriate choice to be discretized and solved. For this 
principle to be used, however, the temporal integration must be removed. This may 
be accomplished through the use of the earlier assumption that the unsteady flow is 
harmonically varying with frequency u. With the temporal behavior of the potential 
represented by assumed mode shapes, the resulting steady and unsteady variational 
principles for the spatial behavior of the flow may then be discretized and solved 
using traditional finite element techniques. 

2.2.2 Steady Flow Variational Principle 

For steady flows, the variational principle is quite simple. Since the pressure, P, 
the potential, $, and the density, R, are all independent of time, the functional II 
becomes 

II steady = JL P dx dy + y> ds (2.60) 

To extremize this equation, the first variation is set to zero, resulting in 


^ri s teady — j j-£ ^ ^ ^ j 

J) Q 6$ ds 


= — j J RV<£> ■ dx dy + j 

f Q <5$ ds — 0 
r 

(2.61) 

Using the divergence theorem, this becomes 



<m ste ady = J t V ' (^V$)] 6$ dx dy + j> 

(q-R^\ 6$ds 

(2.62) 


The Euler- Lagrange equation of this variational principle is the steady conserva- 
tion of mass, Eq. (2.8). The natural boundary condition is the Neumann condition 
Rd$/dn = Q. Finally, it should be noted that the steady conservation of mass is 
nonlinear in the steady potential, $ due to the dependence of R on $. Consequently, 
nonlinear finite element techniques will be required. 

2.2.3 Unsteady Flow Variational Principle With Deforming 
Grid and Vorticity 

The variational principle for the unsteady flow is considerably more complicated than 
its steady flow counterpart. Part of this complexity is due to the fact that the lin- 
earized potential equation is more complicated than the steady flow potential equa- 
tion, as was shown earlier. In addition, we will see that an extension to the variational 
principle is useful for the analysis of blade motion problems. 

Up to this point, the theoretical development in this chapter has been mainly 
concerned with the forced response (or gust) problem, i.e., how to include the effect 
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of unsteady vorticity within a linearized potential framework. There is an additional 
complication inherent in the flutter (blade vibration) problem, however. Previous 
methods for solving the potential equation have used grids that are fixed in space. 
While this is numerically convenient, the airfoil boundary conditions must be applied 
at the mean airfoil location instead of the instantaneous location. As a result, extrap- 
olation terms must be added to the unsteady airfoil boundary conditions to transfer 
the boundary conditions from the instantaneous location of the airfoil to its mean 
position. These terms contain steady velocity gradients that are difficult to evaluate 
numerically. Consequently, the accuracy of numerical calculations using fixed grids 
is limited. 

One way to eliminate these extrapolation terms is to use a computational grid 
that continuously deforms with the airfoil motion. This procedure has been shown to 
be an effective technique for improving the accuracy of linearized analyses [30, 31, 35]. 

In this report, it is assumed that the grid motion is a small harmonic perturbation 
about the mean grid location. We introduce computational coordinates so 

that the grid deforms in the physical coordinate system ( x,y,t ), but appears station- 
ary in the computational coordinate system. Said another way, the computational 
coordinates are “attached” to the deforming grid. To illustrate, Figure 2.3 shows 
contours of the computational coordinate system for deformed and undeformed grids 
for a cascade of airfoils pitching about their midchords. Note that the computational 
coordinates do indeed “deform” with the blade motion. Mathematically, the two 
coordinate systems are related by the perturbation series 


x(Z,ri,T) = £ + f{Cri)e 3 “ T 

(2.63) 

y((,v,T) = v + g(Z,vy wT 

(2.64) 

CS. 

-Cj 

s 

II 

(2.65) 


Note that to zeroth order, the physical and computational coordinates are the same. 
The terms / and g are first-order perturbations that map the computational coordi- 
nate system to the moving physical system. The flow decomposition described earlier 
may be expressed in a similar fashion, i.e., 


KxiV’t) = 


+ 

Pit, l)^ r 

II 

Sh 

<> 

V $(£, g) 

+ 

v (£> 7 ?) eJU,T 

p(x,y,t ) = 

P(Lri) 

+ 

p{t,v) eJWT 

y N) ~~ 

' v- ' 

' v- " 


' — — ' 

total flow 

nonlinear 

mean 


small harmonic 
disturbance 


( 2 . 66 ) 


To first order, the gradient operator, V, and time derivative operator, d/dt , op- 
erators may be expressed in terms of the computational coordinate operators V and 
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Figure 2.3: Contours of computational coordinates (£, 77 ) for undeformed (top) and 
deformed (bottom) grids for a cascade of airfoils. Airfoils are pitching about their 
midchords with an interblade phase angle, a, of 180°. 
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djdr as follows: 


_ ( d/dx \ _ [ 1 - ft -g t 


d/dy 


~fv 1 - 9v 


dm 

d/drj 


= [J]V' 


(2.67) 


d_ _ d_ _ df , 

dt dr dr 


( 2 . 68 ) 


where f is the vector of grid motion perturbation functions, (/, g) T . 

Now that the coordinate systems have been described, the unsteady variational 
principle may be expressed in the computational coordinate system. By assuming 
a harmonic time dependence, the temporal integration in the functional II may be 
carried out. In the computational coordinate system, the unsteady functional is 

Ilunsteady = ^ ^ J J^P |[J] _1 | dr ) Q4> [J] _1 ds dr (2.69) 

The expression for p in the computational coordinate system may be derived by 
substituting the coordinate transformation operators from Eqs. (2.67) and (2.68) into 
the unsteady nonlinear Bernoulli Equation, i.e., 


P — Pt 1 


7- 1 [1 




(2.70) 


Also, to first order, the inverse of the transformation matrix [J] may be written as 


prl— l _ 1 + fi 9t 

LJ ” fv 1 + 9n 


(2.71) 


To obtain an expression for the small disturbance unsteady flow, the functional in 
Eq. (2.69) is expanded in powers of (f>, /, and g retaining up to quadratic terms. 
The first-order terms in the functional will result in the steady flow Euler-Lagrange 
equation. Because the steady flow Euler-Lagrange equation is satisfied by the mean 
flow potential <&, these terms will not contribute to the unsteady flow functional. 
Therefore, only the second-order terms will contribute. Hence, the functional may be 
written as 

II linear = J JJ j^ R [-V^V* + (V'$ T V> + 4 > T ) “] d£ dp dr 

~lLJL R { [ V ' $T ^ + V ' • f (V'$ T V> + </>r) - f T • VV 

~h - f - • v ' $ ) ( v '* Tv '^ + ^)] } ^ dr > dr ( 2 - 72 ) 

Note that the surface integral terms have been omitted for the time being. Next, the 
time integral needs to be evaluated. As was stated earlier, it is assumed that the 
unsteady flow variables (and the blade motion, if any) vary harmonically. Hence, the 
unsteady flow may be represented as a Fourier series in time. These Fourier modes 
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may be analyzed individually and summed, since the governing equation is linear. 
Therefore, it may be assumed that, for example, 


<Kt, V, t) -> Re [<£(£, rjy wt \ = l - [^(£, T})e° wt + «^(£, r,)e J ' wt ] (2.73) 


where <^(£,77) is now the complex amplitude of the unsteady potential, and <^>(£,7 /) 
is its complex conjugate. The other unsteady variables are represented in a similar 
fashion. Substitution of this assumption into the unsteady flow functional, Eq. (2.72), 
results in 


ff linear = \ J R + ~ [V'/VW$ T V> 


+ju> ( V'/ - <j)V§ T + w 2 #] } d£ dr, 


- IJ R { V'$ T [J] + V' ■ f (v'$ T W - - jwf • V> 


1 

c 2 


^ V'$ T [J] V'$ - joj f • ) ( V'$ T W - icu 


d£ dr, 


+ complex conjugate terms (2-74) 

where [J] = [J] T [J] — [I]. Note that none of the boundary conditions have been 
discussed as yet. Additional terms will be added to the variational principle in the 
next section to specify the boundary conditions. 

At this point, it should be noted that no vorticity has been included in the un- 
steady functional. Although the grid motion and unsteady vorticity problems are 
usually analyzed separately in practice, the vortical terms will be added here for 
completeness. 

The divergence theorem says that for any vector V, 



V • V dA- 


i 


V • n ds = 0 


(2.75) 


Consider the case where V = <t>R\ R . The divergence theorem may be written as 


j J V' • Rv r 4> d£ dr, — Rv r • ruj) ds = 0 (2.76) 

Note that because the rotational velocity and grid motion terms in the functional are 
both second order, there is no coupling between the two, i.e., there will be no terms 
in the functional that depend on both the rotational velocity and the grid motion. 
Consequently, Eq. (2.76) may be written in the computational coordinate system. 

We wish to add terms to the unsteady functional II so that the Euler-Lagrange 
equation will contain the appropriate rotational velocity terms appearing in Eq. (2.21). 
Because the terms in Eq. (2.76) sum to zero, Eq. (2.76) and its associated complex 
conjugate expression may be added to the functional given in Eq. (2.74). Leaving 
aside the boundary terms for the moment, taking the variation of Eq. (2.74) includ- 
ing the rotational velocity terms and applying the divergence theorem results in the 
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Euler-Lagrange equation for the unsteady flow. This linearized potential equation 
may be written as 


V' • RV(f> - V' 


L. (v'$ T v'^ + jut) v'$| - A (juV'* T w - l,v) 


= -V' • [i? ([J]V'$ + V' • fV'$)] + jut ■ VR 

r R a 


+v'- 


~ (^V'$ r [j] V'$V'$ - juf • V 


R_ 

+ Q2 


+uH ■ V'$ 


-V'-Rv 


R 


(2.77) 


Note that if there is no grid motion, Eq. (2.77) may be rearranged into the form of 
Eq. (2.21), the original equation to be solved. 

In summary, then, Bateman’s variational principle may be extended to include 
both the flutter and forced response problems typically encountered in turbomachin- 
ery. The next two sections will describe the boundary conditions required for cascades, 
and Chapter 3 will discuss how Eq. (2.77) and the boundary conditions are discretized 
and solved numerically. 


2.3 Near-Field Boundary Conditions 

Now that the governing equations have been developed, appropriate boundary condi- 
tions need to be imposed. Figure 2.4 shows the locations of the four types of boundary 
conditions for cascades. In this report, the periodic, airfoil surface, and wake con- 
ditions will be referred to as “near-field” boundary conditions and will be discussed 
in this section. The remaining “far-field” boundary conditions are somewhat more 
complicated and will be discussed in the next section. 

2.3.1 Periodic Boundary Condition 
Steady Flow 

In this analysis, all of the airfoils in a given blade row are assumed to be identical. 
Hence, the steady flow over each blade is the same. The periodic boundary condition 
is that the difference in the steady potential between two adjacent periodic boundaries 
is a constant, i.e., 

$(£, n + kG) = $(£, v) + kGV- oo (2.78) 

where G is the blade-to-blade gap, V _oo is the specified upstream velocity in the y- 
direction, and k is the blade number, where the reference blade number is zero. This 
boundary condition is applied on the periodic surfaces shown in Fig. 2.4. 
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Unsteady Flow 

Clearly, computing the unsteady flow around every airfoil in a blade row would be a 
formidable computational task. Hence, we wish to reduce the size of the numerical 
calculation. Lane [49] showed that any vibration of the blades of a cascade may 
be decomposed into a sum of traveling waves with a fixed interblade phase angle. 
Accordingly, in this analysis, it is assumed that any incoming disturbance or blade 
motion may be decomposed into a sum of Fourier modes, i.e., the unsteady flow 
is periodic in the circumferential direction. Since the governing equation is linear, 
the flow solutions for each of the individual modes may be superposed to form the 
complete flow solution. 

Hence, the flow around any blade may be computed by simply using the interblade 
phase angle (or circumferential wave number) to account for the distance between the 
blade to be analyzed and a reference blade. The periodic boundary condition may 
then be expressed as 

<Kt,V + *G) = (2.79) 

where a is the interblade phase angle. 

The periodic boundary condition is not included explicitly in the variational prin- 
ciple. Its application will be discussed in Chapter 3. 

2.3.2 Airfoil Surface Boundary Condition 

Steady Flow 

On the surface of the blade, no flow must pass through the blade. For the steady 
flow analysis, there is no blade motion or incoming disturbance, so this condition is 
simply 

*r = ° < 2 - 80 > 

where n is measured normal to the airfoil surface. Note that this condition is the 
natural boundary condition from the steady flow variational principle [Eq. (2.62)]. 
Since it is a natural boundary condition, no boundary condition need actually be 
imposed here to obtain the correct solution. 

Unsteady Flow 

For unsteady flows, the airfoil surface boundary condition is also the natural boundary 
condition of the variational principle. The motion of the blades induces an upwash 
on the airfoil surface (which is the source of the unsteadiness for flutter problems). 
If v R is nonzero on the blade surface (such as in the case of the original Goldstein 
formulation) then the upwash due to the rotational velocity must also be accounted 
for. The natural boundary condition from the variational principle is obtained from 
the first variation of the functional n that includes the rotational velocity terms, i.e., 


^linear = J J^{- ■ •} H d£ df] 
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So the boundary condition is 

d<f) 

dn 


Y ) S<f>ds = ° 

(2.81) 

R -n 

(2.82) 


The first term on the right hand side of Eq. (2.82) is the upwash due to the translation 
of the airfoil. The second term is an additional upwash required to counter a down- 
wash produced by shearing motion of the grid and resulting shearing of the steady 
potential field in the vicinity of the airfoil surface. Note that there is no upwash due 
to local rotation of the airfoil, nor are there extrapolation terms due to the difference 
between the mean and instantaneous position of the blade surface. These are not 
required because the mean potential field moves with the grid due to the grid motion 
assumption imposed earlier. The last term on the right hand side of Eq. (2.82) is the 
upwash associated with a nonzero rotational velocity on the blade surface. 


2.3.3 Wake Boundary Condition 

Steady Flow 

The wake boundary conditions have characteristics of both the periodic and airfoil 
surface boundary conditions described earlier. The first boundary condition is that 
the wake is considered to be an impermeable surface, so the steady no throughflow 
condition applies, Eq. (2.80). Second, there is a jump in potential at the trailing edge 
which is equal to the circulation around the blade. If the jump in potential across 
the wake is equal to the jump in potential at the trailing edge, the Kutta condition 
is automatically satisfied. The boundary condition is then 

[$]wake = ff$]TE (2.83) 

where |$] is the jump in the steady potential. An equivalent boundary condition is 
that the pressure is continuous across the wake. Hence, the steady pressure continuity 
condition is 

[P] = 0 (2.84) 

where [PJ is the steady pressure jump across the wake. For the wake to coincide with 
the grid boundary, an additional equation will be required to enforce the pressure 
continuity condition. 

Unsteady Flow 

For unsteady flows, the wake oscillates harmonically about its mean position with an 
unknown displacement, r. Thus, as in the steady case, there is an auxiliary equation 
to enforce the unsteady pressure continuity. In one dimension, the unsteady pressure 
is linearly related to the unsteady potential, i.e. , 

(2 - 85) 
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The pressure continuity equation may then be written as 

( 3$ d \ 

M = -r (*, + TsT J w = o 


(2.86) 


where [. . .J represents the jump across the wake. Specifically, the jump in potential 
is 

l<t>l = e ~ i£r - 4? (2.87) 

where <f) u is the potential on the upper wake surface of the computational domain (see 
Fig. 2.4), and 4> e is the potential on the lower wake surface. 

The other boundary condition applied on the wake is that conservation of mass 
must be satisfied. In particular, the upwash due to the wake displacement must be 
accounted for. This additional upwash may be written as 


d(f) 3$ dr 

d^ = -dITs +] “ r 


(2.88) 


where d§/ds is the mean flow tangential velocity along the wake. The first term is 
a “ramp” effect due to a change in wake displacement along the wake. The second 
term is the upwash due to the translation of the wake. However, these effects are in 
addition to the grid deformation effects given by Eq. (2.82). Therefore, the complete 
upwash on the wake surface is 


8 $ 

dn 


= juf ■ n — [J]V4> • n + v R 


3$ dr 

n+ 


(2.89) 


Consequently, a term must be added to the variational principle to account for this 
additional upwash due to the motion of the wake. Taking the first variation of the 
modified variational principle and setting to zero gives 

Senear = [Eq. (2.81)] + jf -R + JUT^ - 8^) ds = 0 (2.90) 


2.4 Far-Field Boundary Conditions 

The boundary conditions at the far field are of a different nature than the near- 
field conditions, particularly for unsteady flows. Specifically, the far-field boundary 
conditions isolate the computational domain from influences outside the blade row 
being examined. Unsteady waves from the airfoil must pass out of the computational 
domain unreflected, so that the reflection of these waves does not corrupt the unsteady 
solution. For this reason, it is crucial to apply highly nonreflective far-field boundary 
conditions. 


2.4.1 Steady Flow 

First, the steady flow boundary conditions must be enforced at the far field of the 
computational domain, as shown in Fig. 2.4. The steady far-field boundary conditions 
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are specified using four parameters: the total pressure, pr, the total density, px, 
the inflow circumferential velocity, Vi-oo, and the mass flux through the cascade, 
Qoo — RooUoo- For a fluid with ratio of specific heats 7 , these four quantities may be 
determined once the inflow Mach number and flow angle are known. 

Upstream, the inflow circumferential velocity is used to specify the steady poten- 
tial, $, at the far field, i.e., 

<%) = U_oo y (2.91) 

Note that since this is a Dirichlet boundary condition, the variation of the steady 
potential, <5$>, is zero. Hence, the boundary integral in the first variation of the 
steady variational principle, Eq. (2.62), is zero. 

Downstream, the steady mass flux is is specified, i.e., 

(2.92) 

In this case, a Neumann boundary condition is applied, which is the natural bound- 
ary condition of the steady variational principle. Note that because the governing 
equation is the steady conservation of mass, specifying the mass flux at the exit of 
the domain also fixes the mass flow at the cascade inlet. 

2.4.2 Unsteady Flow 

There are two main approaches for imposing nonreflecting boundary conditions at 
the far field. The first approach, developed by Verdon et a] [50], matched an an- 
alytical description of the far-field behavior to the numerically calculated near-field 
solution. This approach was also applied to the linearized Euler equations by Hall 
and Crawley [29]. The advantage of this method is that the far-field solutions may 
be obtained with very little computational expense, since the far-field behavior may 
be described analytically. There are two main disadvantages to this approach, how- 
ever. First, it does not lend itself well to three-dimensional problems. The annular 
geometry makes analytical solutions extremely difficult to obtain, except for special 
cases. Second, the assumption that the circumferential eigenmodes of the solution 
are Fourier modes is only exact if the computational grid is uniformly spaced in the 
circumferential direction in the far field. In addition, the method does not account 
for the truncation error associated with the field discretization scheme. 

To alleviate these and other problems, Hall [24] and Hall, Lorence, and Clark [51] 
developed a numerically exact method for constructing nonreflecting far-field bound- 
ary conditions. In this approach, the eigenmodes of the unsteady fluid motion in the 
far field are obtained using the discretized governing equations. These eigenmodes 
are then used to construct nonreflecting boundary conditions. Although the numeri- 
cally exact behavior is obtained using this method, it is much more computationally 
expensive than the analytical approach. 

Initially, the numerically exact method was implemented in the present method. 
Unfortunately, the sensitivity analysis (to be described in Chapter 4) of this method 
becomes extremely computationally expensive due to the need to compute the sen- 
sitivities of the eigenvectors. To reduce the computational expense of the present 


Rd i = «■ 
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method, the analytical approach was chosen to demonstrate the present method. For 
the cases to be considered here, the error associated with this approach is actually 
quite small. In the future, particularly for three-dimensional applications, it would be 
useful to find a computationally efficient procedure to perform the sensitivity analysis 
on the exact approach. 

In this section, the method for computing the analytical description of the un- 
steady potential in the far field is described. Only the homogeneous (i.e., the right- 
hand side of Eq. (2.77) is zero) problem will be described here. It should be noted 
that there is never any grid motion in the far field. If there is an incoming acoustic or 
vortical gust, the problem is inhomogeneous, and the particular solution is computed 
numerically. The computation of the particular solution will be described in Chapter 
3. 

Continuous Far-Field Potential 

In the upstream far field, the unsteady potential is continuous in the circumferential 
direction. Downstream, the jump in potential associated with the wake requires an 
additional discontinuous part of the unsteady potential. The continuous part of the 
potential is derived in the same fashion upstream and downstream, and is described 
here. 

In the far field, it is assumed that the steady flow is uniform and the unsteady 
flow is a small harmonic perturbation about the steady flow. Under these conditions, 
the unsteady flow is governed by the convective wave equation, i.e., 

^ - C 2 VVc = o (2.93) 

where the total derivative is 

D_ _ d_ d_ d_ 

Dt dt + dx dy 

Here <f> c is the continuous unsteady potential, U and V are the steady flow velocities in 
the x- and y-directions, and C is the steady flow speed of sound. Next, it is assumed 
that the unsteady potential in the far field may be decomposed into a sum of Fourier 
modes, so that 


4>c{x,y,t) = Y dm exp [jut + ja m x + j/3 m y\ (2.94) 

m=—oo 

where d m are the Fourier coefficients, a m are the as yet undetermined axial wave 
numbers, and f3 m = (a + 27rm)/(7 are the specified circumferential wave numbers. 
The axial wave numbers may be obtained by substituting Eq. (2.94) into Eq. (2.93), 
which gives 

U(u + p m V) ± CyJ(3l(U* + U 2 - C 2 ) + 2f3 m uV + u, 2 

a " m ~~ Q2 _ JJ2 


(2.95) 
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These wavenumbers correspond to upstream going and downstream going pressure 
waves. The direction of propagation may be determined as follows. If the a m pair 
are complex, the direction may be determined on physical grounds by the direction of 
exponential decay. For example, at the upstream far field, the wave that decays in the 
positive x-direction must be a downstream going mode, since if it were an upstream 
going mode, that would imply that the wave generated in the computational domain 
would eventually have infinite amplitude, which is not physically possible. Waves 
that decay exponentially are referred to as “subresonant”. 

If the a m pair are purely real, the waves propagate unattenuated. These are 
referred to as “superresonant” modes. Their direction of propagation may be deter- 
mined by examining the group velocity [51], which is defined as 



(2.96) 


The group velocity is the speed at which a packet of waves of a given wave number 
moves. If the group velocity is positive, there is a net flux of energy in the positive 
x-direction. Hence, the wave propagates in the positive x-direction. 

Finally, the coefficients d m may be determined from a Fourier transform of the 
unsteady potential in the far field, i.e. , 

1 f G 

dm ~ G Jo dy (2.97) 

where x# is the x-coordinate at the far-field boundary being examined. Once all 
of the Fourier coefficients and wave numbers have been obtained, only the outgoing 
pressure waves are retained. All incoming modes are set to zero, since they come 
from outside of the computational domain. The outgoing modes are then summed 
using Eq. (2.94) to determine the new continuous unsteady potential, <f> c . 


Discontinuous Far-Field Potential 

Next, an expression for the discontinuous part of the downstream far-field potential 
due to the wake is required. The wake is modeled here as a vorticity wave that 
convects with the free stream. The geometry for this analysis is shown in Fig. 2.5. 
Here the x- and y - axes form the usual coordinate system, shifted so that the origin 
is at the downstream far-field boundary. The x r - and j/ r -axes are rotated through 
the angle 9 to align with the wake boundary. The wake boundary is aligned with the 
steady free stream velocity, V r . The conversion between the two coordinate systems 
is given by 

x T = x cos 9 + y sin 9 
y r = —x sin 9 + y cos 9 


and 


x — x T cos 9 — y T sin 9 
y — x T sin 9 -fi y T cos 9 
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Figure 2.5: Coordinate systems for downstream far-field analysis. 


Because there is no pressure associated with the wake vorticity, the discontinuous 
part of the potential satisfies Laplace’s equation, i.e., 

VVd = 0 (2.98) 


where (f> d is the discontinuous part of the unsteady potential. 

It will be clearer to derive cf>d in the coordinate system aligned with the flow. The 
vorticity convects with the free stream, so the potential may be assumed to be of the 
form 

<f>d{x r ,y T ) = (j) Ttd {y T )exp(ja T x r ) (2.99) 

where a r = —u)/V r . Substituting this expression into Laplace’s equation gives 


d 2 <f> r ,d{yr) _2j. f_. \ n, 

T 1 a T <p r , d (y r ) = 0 

dy 2 r 

The solution to this ordinary differential equation has the form 

<f>r,d{yr) = Vexp(a T y r ) + £ exp (-a r y r ) 


(2.100) 


( 2 . 101 ) 


There are two boundary conditions required to obtain the coefficients V and £. The 
first is that there is a jump in potential across the wake, which is defined as 


l<f> T j - <f>r,d{ G r) exp(- j&Gv ) - ^ r> «*(0) (2.102) 
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where fi T may be determined from the interblade phase angle, <7, and the coordinate 
transformation 

a + cjGsin 9IV T 

* - Geos# ( 2 ' 103 ) 

The second boundary condition is that the wake motion is periodic 


<t>r,d(0) = 4>T,d{Gr) exp(—j/3 r G r ) 


(2.104) 


After using these two boundary conditions, some algebra results in the expression for 

$r,d 


<t>T,d(yr) = 


exp (a T y r ) 


+ 


exp (~a T y r ) 


(2.105) 


[1 - exp(o; r 6 ? r - j(3 r Gr) 1 - exp(-a r G> -jj3 T G T ) 

The next task is to determine [<^ r |. Using the notation of Figure 2.5, it is clear that 


IM = C exp(-j/3 r G r ) - 4> e - 
Since the behavior of <f>d is known in the x r direction, 

4% = 4> u exp (ja r G sin 9) 

Substituting this expression into Eq. (2.106) and simplifying results in 

{4> t \ = <f> u exp(-jcr) - 

So the complete expression for the discontinuous unsteady potential is 

<f>d{x T ,yr ) = -^Urj exp {ja T x T ) 


(2.106) 


(2.107) 


(2.108) 


x 


exp(a r y r ) 


+ 


exp {-a T y r ) 


1 — exp(o; r G' r — j(3 r G r ) ' 1 — exp(— a r G T — j(3 r G r )\ 

Finally, for this expression to be used in the present method, it must be converted to 
the original coordinate system, so that 

4d(x, y) = M exp [ja T (x cos 9 + y sin 9 )] 


exp[a r ( — x sin 9 + y cos 0)] exp(— a T (— x sin 9 + y cos 9)} 

1 - exp(cx r G r - j/3 r G r ) 1 - exp(-a r G r - j(3 r G T ) 


(2.110) 


This expression is added to the continuous downstream unsteady potential to deter- 
mine the complete solution at the downstream far field. 



Chapter 3 

Numerical Solution Method 


In this chapter, the methods for numerical solution of the governing equations and 
boundary conditions described in Chapter 2 for the steady and small disturbance un- 
steady flow through a two-dimensional cascade of airfoils are presented. Section 3.1 
contains a description of the steady solution procedure. The first step in the steady 
solution process is the choice of a computational grid upon which to discretize the 
governing equations and boundary conditions. After the equations governing the gen- 
eration of the computational grid have been discussed, the finite element discretization 
of the steady flow variational principle and the associated boundary conditions will 
be presented. To complete the steady flow analysis, the matrix assembly and solution 
of the discretized equations and boundary conditions will be described. In Section 
3.2, the unsteady solution procedure will be developed, beginning with the numerical 
calculation of the drift function and rotational velocity, as well as the generation of 
the grid deformation. Next, the unsteady finite element discretization of the field 
equations and boundary conditions is presented, with particular emphasis on the nu- 
merical discretization of the far-field boundary conditions. Finally, the assembly and 
solution procedure for the discretized unsteady flow equations is described. 


3.1 Steady Flow Solution Procedure 

3.1.1 Grid Generation 

In any computational fluid dynamic analysis, the choice of a computational grid 
upon which to discretize the governing equations is crucial. Poor computational 
grids may reduce the accuracy of the discretization scheme or, in some cases, reduce 
the stability of the solution procedure. For cascade flows, the computational grids 
typically used are referred to as C, O, and H grids (or some combination thereof). 
The grid topology for each of these grids looks much like the letter used to identify 
them. In this investigation, H grids will be used. H grids, unlike their C and 
O counterparts, provide good resolution throughout the computational domain, not 
just near the airfoil. This is important because the accurate resolution of acoustic 
and vortical waves is critical to evaluating the aeroacoustic performance of a blade 
row. Another important reason for using an H grid is that the matrix containing 
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y 



-> X 


Figure 3.1: Typical computational grid for a fan exit guide vane. 

the discretized flow equations will be block-tridiagonal. In addition to reducing the 
memory requirements for the method, block-tridiagonal matrix solvers are, in general, 
considerably faster than a general (i.e., fully populated) matrix solution. A typical 
H grid is shown in Figure 3.1. 

There is an additional requirement of the computational grid. Recall from Chapter 
2 that for aeroacoustic calculations, knowledge of the drift and stream functions is re- 
quired. To facilitate the computation of these functions, one possible approach would 
be to solve for the steady flow on some initial computational grid, and afterwards cal- 
culate the streamlines based on the computed steady flow. Hall and Verdon [28] used 
this approach, finding the stagnation point and then calculating the stream function 
using a Runge-Kutta algorithm. Unfortunately, this procedure increases the overall 
truncation error of the solution procedure and does not lend itself well to a sensitivity 
analysis. A method that is better suited to a sensitivity analysis is to compute the 
grid nodes as part of the steady solution procedure. The final steady solution will 
then be on a computational grid that follows the streamlines. The latter procedure 
will be implemented here. 

Thompson [52] has developed an elliptic grid generation technique that is well- 
suited for cascades. The ( x,y ) position of the grid nodes are defined by the partial 
differential equations 

V 2 E'= 


V 


(3.1) 
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V 2 H = Q 


(3.2) 


where V and Q are functions which are used to control the grid spacing. Lines of 
constant E and H define the grid lines. 

Equations (3.1) and (3.2) may be inverted to obtain partial differential equations 
for the unknown grid node locations in terms of the known computational coordinates 
E and H , i.e., 

_ d 2 x — d 2 x _d 2 x - 2 -da; 

b ~ ar7 + 7TTTT + 8 Q~ = 0 


<9E : 


dEdH ' dH 2 


'dH 


and 


where 


_d 2 y — d 2 y _d 2 y -^2 ~ dy 

“as 5 " 2l3 asdH + + 6 q bh = 0 


(3.3) 

(3.4) 


a = 


' dx' 

dH 

\ / 

dx 





dxdx dy dy 
^ = ~dEdH + dEdH 

— dx dy dx dy 
= 9E~dH ~ dHdE 


The next task is to choose the grid spacing functions V and Q so that the grid 
follows the flow streamlines. These functions may be obtained through the fluid dy- 
namic definitions of the velocity potential and stream function. Consider an inviscid, 
irrotational, compressible flow. In terms of the stream function, 4', and the density, 
R, the steady conservation of mass may be written as 


d_ 

dx \ R dx 


d_ 

dy V# dy t 


= 0 


(3.5) 


which, after rearranging, may be written in the form 

T ldVdR 1 dtf dR 

V 2 ']/ = 1 

R dx dx R dy dy 


(3.6) 


Finally, writing the right hand side of Eq. (3.6) in terms of the steady velocity po- 
tential, $, gives 

V 2 ^ = |V$ x Vi?| (3.7) 

Note the similarity between Eq. (3.7) and Eq. (3.2). If we choose V = 0 and 
Q = |V$ x Vi?|, then lines of constant H will correspond to streamlines, provided 
the boundary conditions around the computational domain are consistent with the 
definition of the stream function. Furthermore, note that if the flow is incompressible, 
| V$ x Vi?J =0, and no grid spacing functions are necessary in the grid generation 
equations (i.e., V — Q = 0). Regardless of the compressibility of the flow, it should 
be noted that lines of constant E do not correspond to contours of the drift function, 
A, nor is this a requirement of the computational grid. The drift function may be 
calculated after the steady flow solution has been obtained. 

Figure 3.2 shows a typical computational grid in the (E, T) coordinate system. 
If the grid generation equations, Eqs. (3.3)-(3.4), were solved without the influence 
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Figure 3.2: Typical computational grid in coordinates. The computational 

node numbering convention is also illustrated. 


of the grid spacing functions (i.e., V = Q = 0), an example of the resulting compu- 
tational mesh is shown in Figure 3.1. Note that the area of the computational cells 
shown in Fig. 3.1 varies smoothly throughout the computational domain. The lack 
of large area changes from cell to cell is characteristic of elliptic grid generation al- 
gorithms, and is desirable because smoothly varying grids do not increase the overall 
truncation error of the numerical solution procedure. 

The overall computational grid has I nodes in the axial direction and J nodes 
in the circumferential direction. The node numbering convention is illustrated in 
Figure 3.2. Any node on the computational grid may be referred to by its node 
number (i, j). The upstream far-held boundary is the first axial grid line, so i = 1 on 
this boundary. Similarly, i = / on the downstream far-held boundary. The j-nodes 
are numbered from the lower boundary of the computational domain ( j = 1) to the 
upper boundary of the domain (j = J)- 

In addition to the grid generation equations described above, boundary conditions 
must be imposed around the blade passage. Specifically, the distribution of grid nodes 
along the boundaries are specihed. These locations are specihed as a fraction of arc 
length on each boundary. There is no restriction on the grid distribution, although to 
minimize truncation error, it is best to increase the grid resolution near large steady 
how gradients. A typical boundary distribution is shown in Figure 3.3. The arc 
length fraction for each boundary node i along the periodic, airfoil surface, and wake 
boundaries is denoted by Fi. The value is the same for both the upper and lower 
surfaces. The arc length fraction for each node j along the far-held boundaries is 
denoted by Fj. The values of Fi and Fj at the endpoints of each boundary type are 
also shown in Figure 3.3. The values of Fi and Fj at the intermediate nodes represents 
the fraction of the arc length between the endpoints at which the node is located. 

The grid generation equations and associated boundary conditions are discretized 
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Fi =0 


Fi = 1 


Fi = 0 



x 


Figure 3.3: Typical distribution of boundary grid points. 

using centered finite difference operators for nonuniform spacing in E and tH. The 
resulting set of nonlinear equations is of the form 

M($,x;Z) = 0 (3.8) 

where M is a vector of nonlinear functions, $ is the vector containing the discrete 
approximation of the nominal steady potential $ at all the computational nodes, x is 
the vector of the location of the computational nodes, and Z is a vector containing a 
set of parameters that define the airfoil shape and cascade geometry. The semicolon 
indicates that $ and x are the dependent variables to be computed during the solution 
procedure. The parameters contained in the Z vector are independent variables during 
the solution procedure, although they may change in the overall design process. The 
details of the parameters in the Z vector will be discussed in the next section. 


3.1.2 Airfoil Definition and Spline Notation 

Now that we have defined the computational grid, we consider how to represent the 
airfoil shape and cascade geometry using the vector Z. Generally, airfoil shapes are 
defined using two distinct approaches. One approach is to define the airfoil analyti- 
cally using some set of design variables such as thickness and camber. In this case, 
the vector Z may be considered to contain the magnitudes of each of these defining 
variables, in addition to the cascade parameters. A second, more general approach is 
to define the airfoil shape using a set of points on the airfoil surface. In this case, the 
vector Z will contain the x- and y-coordinates of the defining points as well as the 
cascade parameters. It should be noted that this second approach is more general 
because any analytical airfoil definition may be discretized to obtain a set of defining 
points. Furthermore, it may be argued that each defining point is actually a design 
variable of the airfoil shape. Hence, for the most part, the discussion of the airfoil 
shape in this report will assume that the blade has been defined by a set of points on 
the surface. 
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Figure 3.4: Detail of airfoil spline definition. Arrows indicate direction of airfoil 
spline. 


The cascade parameters are design variables that are independent of the airfoil 
shape. For example, the blade-to-blade gap, G , and the stagger angle, 0, are both 
parameters that would be included in the vector Z. In general, any parameter that 
defines the cascade and is considered to be a design variable (i.e., the parameter 
is allowed to change during the design process) is included in Z. This may include 
parameters that define the computational grid, such as the fractional arc length arrays 
F, and Fj defined in the previous section, and parameters that define the steady flow 
conditions, such as the inflow Mach number. 

At this point, it is useful to explain the relationship between the vector containing 
the airfoil shape and cascade parameters, Z, and the vector containing the nodes of 
the computational grid, x. Consider a blade definition consisting of a number of 
defining points as shown in Figure 3.4. Each point defining the shape of the blade 
has a coordinate that will be denoted (X, Y). In the present analysis, it is assumed 
that the points begin at the trailing edge on the upper surface of the airfoil and end 
at the trailing edge on the lower airfoil surface. Note that in general these points are 
not points defined by the computational grid. The coordinates of the grid nodes (i.e., 
the vector x) will be denoted using lower case letters, i.e., ( x,y ), to distinguish the 
nodes of the computational grid from the points defining the airfoil shape. 

Figure 3.4 shows a schematic of an airfoil defined by some finite number of points, 
N. The points begin at (Xj, Y\) and end at (Xn, Fn). In addition, the initial estimate 
of the stagnation point is denoted by (Xstag, Fstag)- The arc length, S, of each of 
these points measured from the trailing edge of the upper airfoil surface is then 
calculated at each point based on its X and Y location. Once the arc length has been 
calculated, the X and Y locations are given a “functional” definition using a cubic 
spline. In this way, the X and Y location of the surface of the airfoil is defined to be 
a function of S. 

Although the airfoil shown in Fig. 3.4 has a sharp trailing edge, it should be noted 
that this analysis is not restricted to airfoils with sharp trailing edges. Actual airfoils 
usually have rounded (or blunt) trailing edges. Unfortunately, the wake boundary 
condition is difficult to apply near rounded trailing edges. In such cases, it is useful 
to “cut off” the airfoil near the trailing edge so that there is a space between the first 
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Figure 3.5: Typical computational grid illustrating modifications for airfoils with 
blunt trailing edges. Multiple passages shown for clarity. The highlighted lines define 
the actual computational grid. 

and last point of the airfoil definition. The resulting computational grid around such 
an airfoil is shown in Figure 3.5. Note that there is constant spacing between adjacent 
blade passages from the trailing edge to the downstream far field of the computational 
grid. This space is intended to model the flow deficit due to a viscous wake. Placing 
the grid boundaries on either side of this space is analogous to representing the 
displacement thickness of a viscous boundary layer by a solid surface in an inviscid 
analysis. Although no displacement thickness modification has been applied to the 
airfoil surface, the space in the wake represents an estimate of the actual viscous wake 
region. Inside the space between the grids there is no fluid flow. The wake boundary 
conditions are then applied in the usual fashion, on either side of the space. Previous 
researchers [35, 24] have shown that modeling a blunt trailing edge through a space 
in the computational grid results in a solution that satisfies the Kutta condition, and 
hence no spurious stagnation points arise along the wake boundary. 

Finally, there are four other parameters that define the computational grid. The 
first two are the x-location of the upstream far field, xmin, and the inflow angle fl_oo- 
The other two are the x-location of the downstream far field, xmax> and the flow exit 
angle, fioo- Although the inflow and exit flow angles may change during the course of 
the steady solution procedure, xmin and xmax remain fixed. These four parameters 
are shown in Fig. 3.5, and they complete the specification of the computational grid. 
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3.1.3 Finite Element Discretization 

The next step in the analysis is to calculate the steady flow. As was shown in 
Chapter 2, the governing equation is the steady version of the full potential equation, 
Eq. (2.11). Following the notation in the previous section, the discretized set of steady 
flow equations has the form 

N($,x;Z) = 0 (3.9) 

where N is a vector of nonlinear functions. Note that Eqs. (3.8) and (3.9) are both 
nonlinear in $ and x. Since the objective given earlier is to solve the equations 
simultaneously, Newton iteration will be used to solve the grid and flow equations. 
Consequently, it will be useful to discretize the steady flow equations in a manner 
consistent with Newton iteration. 

From the functional given in Eq. (2.62), it is clear that a perturbation expression 
for the steady pressure is required. To do this, consider the steady Bernoulli equa- 
tion, Eq. (2.10). If the current estimate of the steady pressure is P and the current 
current estimate of the steady potential is $, then to compute the new value of the 
steady pressure, we expand Eq. (2.10) in a Taylor series about the current value of 
P, retaining terms to second order. In terms of the perturbation potential, the 
new estimate of the steady pressure may be written as 

P N E w = P~ PV$ • W - f V$' • V$' + (V$ • V$') 2 + 0 ( $' 3 ) (3.10) 

z zO z 

where Pnew is the new estimate of the steady pressure. In Chapter 2, it was noted 
that Taylor expansions within variational principles must be carried out to second 
order if a first-order result is desired. This is because taking the variation will reduce 
these second-order terms to first order. Since first-order terms are required in a 
perturbation expression for the pressure, the Taylor expansion in Eq. (3.10) is carried 
out to include terms of second order. 

Now, rewriting the steady flow functional given in Eq. (2.62) to second order 
results in 

rigteady = // £ [p ~ PV$ • V$' - |v$' • V<N + ~ (V# • V*')’] # (3.11) 

where the boundary integral has been omitted for the time being. Taking the first 
variation and setting to zero gives 

^n st eady = jj R [- • V<J$' - V$' • 

+ ^-V$ • V$'V$ • V£$ , l d£ drj = 0 (3.12) 

For II s teady to be stationary, 6n ste ady must be zero for all admissible variations in 
Since the area of the domain E is arbitrary, the expression inside the integral 
in Eq. (3.12) must be zero (i.e., it is the Euler-Lagrange equation of the variational 
principle). Hence, discretization of the functional given in Eq. (3.11) will result in a 
set of nonlinear equations of the form given in Eq. (3.9), which is the desired result. 
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The variational principle may be discretized using conventional finite element 
techniques. In the present work, a four node bilinear isoparametric element will be 
used. The values of the perturbation steady potential at the corners of a given element 
n may be interpolated into the interior of the element using an interpolation of the 
form 

m,V)= |NJn*'n (3.13) 

where 

|NJ n = [N u N 2 ,N 3 ,N 4 \ n (3.14) 

is a row vector of interpolation functions and $' n is the vector containing the corner 
values of the perturbation steady potential. Similarly, the gradient operator may be 
represented by the gradient of the interpolation functions, i.e., 


where 


[N']„ = 


V3>' = [N'„]$' n 

dfVy dNji dNz dN± 

d( 5 d( ’ di ’ di 


(3.15) 


(3.16) 


a/V]L dN^ dN± dN A 
dr} 5 dr} 5 dr/ ’ dr} 

Substituting these expressions into the functional given in Eq. (3.11) results in 

R 


n s teady = / jf [p ~ R**[ N'£Vd> - f <[N'£[ N %$' n 


R 


2 C 2 


* # ;[N1Iv*v* T [Ni w *' T 


d( dr} 


(3.17) 


Setting the first variation to zero results in the Newton iteration equation for the 
perturbation of the flow 


+^t[n'£v$v* t [N')„s'„ 


sn, u . dy = JJ^ [-[N'ffv* - [N'fflN'].#'. 

R dx dy — 0 ( 3 . 18 ) 

j 

Rearranging, this may be written as 

a, = / [[N'£ (-[i] + iviv# r ) [N'] n *' n 

-[N'£V$] R d£ dr} = 0 ( 3 . 19 ) 

In finite element notation, then, the steady flow equations for each element may be 
written in the form 

[K] n $'„ - E n = 0 (3.20) 

where the elemental stiffness matrix, [K] n , is 

[K]„ = jj ^ [N')I (-[I) + ±V.V^) [Nl.fi di d v ( 3 . 21 ) 
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and the elemental force vector, E n , is 

En = // [N'£v$i? d( dp (3.22) 

The elemental stiffness matrix and force vector are computed for each cell in the 
computational domain. These elemental values are assembled into a global stiffness 
matrix and force vector so that the steady flow governing equation, Eq. (3.9), may be 
solved using Newton iteration. Also included in the global stiffness matrix and force 
vector are contributions from boundary integral terms resulting from the steady flow 
variational principle. The boundary conditions will be discussed in the next section, 
followed by a description of the global matrix assembly and solution procedure. 


3.1.4 Near-Field Boundary Conditions 

Periodic Boundary Condition 


In Chapter 2, it was stated that the steady periodic boundary condition is that the 
difference in the steady potential between two adjacent periodic boundaries is a con- 
stant. In addition, it is clear from the discussion of the grid generation procedure 
that the periodic surfaces on the upper and lower boundaries of the computational do- 
main must be separated by the specified blade-to-blade gap, G. There is an additional 
requirement of the steady flow solution on the periodic boundaries, however. . 

Earlier in this chapter, we stated that we wish to have the computational grid 
follow the flow streamlines. Because the periodic boundaries are attached to the 
airfoil, the periodic boundaries must be stagnation streamlines of the steady flow. 
Consequently, we wish to formulate boundary conditions to enforce this requirement. 

The stagnation streamline condition may be enforced through the application of a 
“Kutta condition” requirement on the periodic boundaries. This may be accomplished 
by implementing two separate boundary conditions. First, there is a condition that 
is purely a function of the grid locations. In general, the shape of the stagnation 
streamline from the inlet to the leading edge of the airfoil is not known a priori. 
Often the boundary has some curvature, because the stagnation streamline must be 
normal to the surface of the airfoil at the stagnation point, which is not necessarily 
the same direction as the specified flow direction in the upstream far field. Hence, 
the grid locations on the periodic boundaries are not fixed, but are adjusted during 
the solution procedure, so that the same relative spacing along the boundaries is 
maintained. Specifically, on the lower ( j = 1) periodic boundary of the computational 
domain, the boundary condition for the grid is that at grid station i, the x- and y- 
coordinates are the appropriate fractional arc length between the upstream far-field 
boundary and the airfoil surface, so that 


mi = (a+i - F.f [(*, - 

- (Fi — Fi-i) 2 [(x,+i -x. 


^-l ) 2 + (yi - Vi- 1 ) 2 


) 2 + (y.-H - yi? 


= o 


(3.23) 


where the superscript x refers to the x grid equation for the computational node ( i,j ) 
indicated by the subscripts, and F,- refers to the fractional arc length array described 
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in the previous section. This boundary condition is applied from node i = 2 to 
node i = ILE — 1 (see Figure 3.2). At the inlet (i = 1) boundary, the x-location is 
prescribed. The y-location is permitted to change to satisfy Eq. (3.23). 

On the upper (j = J) periodic boundary of the computational domain, the x- 
coordinate of the grid node is the same as the lower boundary, and the y-coordinate 
is one blade-to-blade gap larger than the lower boundary, i.e., 

M*j = x^j - x iA = 0 (3.24) 

and 

M?j = y itJ - y itl - G = 0 (3.25) 

The other boundary condition for closure is purely a function of the steady po- 
tential. In Chapter 2 during the discussion of the steady wake boundary condition, 
it was noted that if the jump in steady potential across the wake is the same as the 
jump in potential at the trailing edge, the Kutta condition is automatically satisfied. 
This condition may be applied in an analogous fashion on the upstream boundary. If 
the jump in potential across the periodic boundary is the same as the jump at the 
stagnation point (i.e., zero), then the periodic boundary is a stagnation streamline of 
the flow. This condition may be expressed as 

= ($ t+ i,j - $<+i,i) - (*.V ~ $.,i) = 0 (3.26) 

The result of these boundary conditions is that the upstream periodic grid boundaries 
are stagnation streamlines of the converged solution. 

Airfoil Surface 

As was noted in Chapter 2, on the airfoil surface, the steady flow satisfies the natural 
Neumann boundary condition, i.e., no through flow. Since it is the natural bound- 
ary condition of the variational principle, no additional expressions are required to 
enforce this condition. The grid equations, however, are another matter. For each 
computational node on the airfoil surface, there are boundary conditions that must 
be satisfied. 

The first set of boundary conditions requires that the computational nodes lie on 
the airfoil surface. For the lower (j = 1) boundary of the computational domain, this 
condition may be expressed by the equations 

M? A =x itl -Xs(si,i) = 0 (3.27) 

and 

Ml=yi,x-Ys{si, i) = 0 (3.28) 

where Xs and Ys are cubic spline evaluations of the airfoil surface coordinate corre- 
sponding to the arc length Si t \ for the computational node at the zth station on the 
j = 1 surface. Note that s, f i is a new variable that is stored only for the airfoil surface 
points (i.e., from i = ILE to i = ITE). This new variable will require an additional 
boundary condition to close the system, which will be described shortly. 
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It should be clear now why the airfoil definition points (A", Y) are referred to 
separately from the grid nodes (x,y). At each iteration of the Newton solver, the 
nodes of the computational grid are permitted to “slide” along the surface specified 
by the airfoil defining points. The coordinates of the points containing the airfoil 
definition are fixed during the steady solution procedure. Hence, these points must 
be considered separately from the grid node locations. 

On the upper (j = J ) surface of the computational domain, the grid nodes must 
also lie on the airfoil surface, so that 

= x itJ - X s (si,j ) = 0 ( 3 . 29 ) 

and 

MIj = Vi,j - YsM - G = 0 ( 3 . 30 ) 

where the latter equation also ensures that the blade-to-blade gap is constant along 
the blade, and s t| j refers to the arc length of the computational node i on the j — J 
surface. 

The other boundary condition for closure is that the arc length that the nodes 
on the airfoil surface “slide” is a specified value of the fractional arc length between 
the stagnation point and the spline endpoint. This condition is applied at all airfoil 
surface nodes except the stagnation point. Recall that earlier in this chapter, the 
arc length fraction F x was defined. Using this definition, the arc length at station i 
and (i — 1) on the lower (j = 1) airfoil surface of the computational domain may be 
expressed as a function of the arc length at the stagnation point and the arc length 
at the trailing edge, so that 


M = — + Ssp — Fi ( Ssp — Si) = 0 (3.31) 


M = — 5,-ip + Ssp — F-i (Ssp — Si) = 0 (3.32) 

where Ssp is the arc length at the stagnation point. Solving both of these equations 
for Ssp and setting them equal to each other results in an expression for 3,-p that is 
only dependent on the specified fraction array F, the arc length at an adjacent node, 
Si-ip, and the arc length at the trailing edge, S\ (see Fig. 3.4). Hence, the boundary 
condition that enforces the fractional arc length distribution on the (j = 1) airfoil 
surface may be written as 


M* J+ a 


l — Fi Fi — Fi- 1 a n 

= 5t ’ 1 ~ Y^fu ” T^f~ 7 ~ 0 


(3.33) 


Note that the J + l designation on M indicates that this equation is solved in addition 
to the grid and flow equations for the J nodes at each i station. Hence, Eq. (3.33) is 
stored in the x grid equation for the “fictional node” J + l. This boundary condition 
is applied from node i = ILE + 1 to i = IT E on the j — 1 surface. 

In a similar fashion, the boundary condition on the j = J airfoil surface may be 
written as „ „ „ 

- Fi- j 


MJ 


,j+i — S *,J 


1 -Fj 

1 — Fi--[ 


Si-1, j — 


1 ~ F-y 


Sn = 0 


(3.34) 
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This equation is stored in the y grid equation for the “fictional node” J + 1 . Like the 
j = 1 surface, this boundary condition is applied from node i = ILE -f 1 to i = ITE 
on the j = J surface. 

These last two boundary conditions essentially state that any change in position 
of the stagnation point will result in a proportional change in the position of the grid 
nodes on the airfoil surface. 


Stagnation Point 

Next, we consider the boundary condition at the stagnation point. This boundary 
condition is extremely important because the description of the unsteady vorticity 
given in Chapter 2 is only valid if the stagnation point is identified accurately. The 
boundary condition is that the steady potential, $, has a minimum value on the 
surface of the blade at the stagnation point. In other words, the derivative of the 
steady potential in the tangential direction (i.e., with respect to the arc length) is zero. 
Discretizing this derivative using a second-order finite difference operator allowing 
variable grid spacing results in 




1 SlLE - SILE+\,J 

SILE+1,1—SILE+1,J |_S/££+l,l “ S ILE 


($ILE+\,J ~ V-ooG - $ILE, l) 


j &ILE+ 1,1 ~ SILE 
SILE ~ s ILE+l,J 


($ILE, 1 ~ $/££+!, l) 


= 0 


(3.35) 


where ILE refers to the grid station i corresponding to the leading edge of the airfoil 
(see Fig. 3.2). The i = ILE grid station is also the stagnation point. Note that the 
potential on the upper ( j = J ) surface of the computational domain must be evaluated 
on the reference airfoil, so the jump in potential between adjacent blades, VL^C?, must 
be subtracted from the value of the upper surface potential, $rLE+i,J- Furthermore, 
note that s only has a single subscript at the leading edge. This is because the arc 
length at the leading edge is the same for both airfoil surfaces. Hence, only one 
additional variable must be stored, which results in there being only one additional 
boundary condition at the leading edge stagnation point. 


Downstream Wake Boundary 

The last near-field boundary conditions are applied at the downstream wake bound- 
ary. Like the upstream periodic condition, the grid line attached to the trailing 
edge must be a stagnation streamline. Hence, the boundary conditions on the wake 
surfaces are nearly identical to those on the periodic surfaces. First, the boundary 
condition on the grid for the lower (j = 1) surface of the computational domain is 


Mli = (F i+ 1 - Fif [(xi - Xi-i ) 2 + (y t - yi- 1 ) 

- (F t - i) 2 \(x i+1 - Xif + (y i+l - yif \ = 0 


(3.36) 
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Also like the periodic boundary, the grid boundary conditions on the upper (j = J) 
surface are 

= x hJ - x iA = 0 (3.37) 

and 

Mi,j = Vi,j ~ Vi,\ - [G — (Fi — Vjv)] = 0 (3.38) 

where the last two terms in Eq. (3.38) account for any space in the downstream wake 
due to a blunt trailing edge, as described earlier in this chapter. 

Finally, similar to the upstream periodic boundary, the boundary condition that 
enforces the Kutta condition may be written as 

= ($i,J ~ $i,i) - = 0 (3.39) 

there are two differences between the wake boundary condition expressed in Eq. (3.39) 
and the periodic boundary condition given in Eq. (3.26). First, the jump in the 
potential is measured at the i and (i — 1) points instead of the i and (* + 1) points. 
This is because the wake boundary conditions are applied from node i = ITE + 1 
to i = I . Hence, so that a single condition may be used at all of these nodes, the 
difference is computed in the “upwind” (i.e., negative i) direction. Upstream, because 
the boundary condition is applied at the i = 1 node, the difference is computed in 
the “downwind” (i.e., positive i) direction. The second difference between the wake 
condition and the periodic condition is that there is a finite jump in potential across 
the wake if there is circulation around the airfoil. 

All of the near-field boundary conditions necessary for the calculation of the steady 
flow have now been presented. The final step in the steady flow solution procedure 
is to examine the steady far-held boundary conditions. 

3.1.5 Far-Field Boundary Conditions 

Upstream Far-Field Boundary 

There are three boundary conditions the grid and steady how must satisfy on the 
upstream far-held boundary. The hrst is that the y-component of velocity must be 
equal to a specihed value. This condition effectively specihes the steady potential 
at the upstream boundary, as shown in Eq. (2.91). Computationally, this may be 
accomplished using the hnite element procedure known as the penalty method [53]. 
Using the penalty method, the boundary condition may be written as 

N hj = [Kl u *, - E, j - k r (i hi - V.„ yiJ ) = 0 (3.40) 

where the matrix [K] and vector E are the stiffness matrix and force vector described 
earlier, and k p is called a penalty number. As k p approaches infinity, the condition 
j = V—oo yi j is exactly enforced. For computational purposes, it is sufficient to 
choose k P to have a large hnite value. 

The other two conditions apply to the location of the grid nodes on the upstream 
far-held boundary. The boundary is a line of constant z-coordinate, while the y- 
coordinate is specihed as a fraction of arc length through the F array described 
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earlier. The ^-coordinate is computed through the equation 

Mlj = - x lfl = 0 (3.41) 

Note that the value of Xij is computed as part of the periodic boundary conditions 
described in the previous section. The y-coordinate is computed in a similar fashion, 
i.e., 

M h = yij - 2/i,i - F i G = 0 (3-42) 

where F } is the value of the far-held arc length fraction array corresponding to node 
j (see Fig. 3.3), and G is the blade-to-blade gap. 

Downstream Far-Field Boundary 

The downstream how boundary condition specihes the mass flux through the compu- 
tational domain. Because the governing equation of the steady how is the conservation 
of mass, specifying the mass flux at the exit hxes the mass flux at the inlet. The cir- 
cumferential grid lines are modeled as streamtubes, so that the how through each 
tube is proportional to its fraction of the gap at the far-held boundary. It should be 
noted that the downstream gap is not necessarily the same as the blade-to-blade gap, 
G, due to the space left in the wake for airfoils with hnite thickness trailing edges. To 
prescribe the mass hux at the boundary, the desired mass flux is added to the “force 
vector” computed from the hnite element procedure, i.e., 

1 G 

N U = - -^R-ocU-oo Q _^y x Ly n ) ( yi ' j+1 - Vl '^ = ° ( 3 - 43 ) 

and 

1 G 

Ni,j+ i = — E/j+i — ^R-ooU {yi,j+i — Vi,j) = 0 (3-44) 

The last term in Eqs. (3.43) and (3.44) represents the mass hux through the stream- 
tube dehned by two adjacent streamlines (i.e., lines of constant j). The grid equations 
at the downstream far-held boundary are analogous to those at the upstream bound- 
ary, i.e., 

Mf j = x/j - x IA = 0 (3.45) 

and 

Mjj = y,,i - vi, i - r,[G - (y, - r„)l = o (3,46) 

where again the blade-to-blade gap, G, is modihed by the space left in the wake due 
to hnite thickness trailing edges. 

3.1.6 Assembly and Solution 

The hnal step in describing the steady how solution method is to explain how the 
equations given in the previous sections are actually assembled and solved. Earlier in 
this chapter, it was noted that the governing equations and boundary conditions for 
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the grid and steady flow could be expressed as two sets of nonlinear functions, M and 
N [see Eqs. (3.8) and (3.9)]. It was also noted that the equations would be solved 
simultaneously using Newton iteration. Schematically, the steady solution procedure 
may be written in the form 


$ \ 

n+1 

H* 

r 5N 

aN 1 
ax 

x ) 


U J 

aM 

aM 


L a$ 

ax J 


-l 


n 



(3.47) 


where n denotes the current estimate of the solution and n + 1 denotes the new esti- 
mate. Literally structuring the equations in this way would result in a sparse matrix 
of extremely high bandwidth. Large bandwidth sparse matrices are not conducive 
to efficient solution. Fortunately, the bandwidth of the matrix may be substantially 
reduced by reordering the equations. 

At each computational node, there are three equations to be solved, for variables 
$, x , and y. Examination of the equations derived in this chapter shows that the 
discretized form of the governing equations and boundary conditions at a grid station 
i are at most dependent on variables at the (i — 1), i, and ( i -f 1) stations. The 
boundary conditions may depend on variables at any j station, however. Hence, the 
equations may be reordered so that the global matrix is block-tridiagonal, where each 
block equation represents a single grid station i. 

For a computational grid with I nodes in the axial direction and J nodes in the 
circumferential direction, the matrix to be solved at every iteration has the form 


' Bj C a 
A2 B2 C2 


A, B, C, 


(3.48) 


A/_i B/_! C/_i 
A/ B / 


where each block row corresponds to an axial grid line. For each i block equation, A,- 
contains terms that multiply variables at the ( i — l)-station, B, contains terms that 
multiply variables at the ^-station, and C, contains terms that multiply variables at 
the ( i + l)-station, so that 


A, 


1— 1 
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-yTl-t-1 

X i-1 


— x ;• 


+ B, 


/ $? +i - 


I -yrP 

\ "S 


+ c,- 


$! i + 1 _ a> 

^1+1 v 
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i+1 


-y-71+1 

A + l 


Y^ 1 

A t+i 



(3.49) 

where [^" +1 — $”,x[ l+1 — x"J T is the change in the solution vector at the i-station, 
and LMi,NjJ r is the residual at the f-station. 

The sub-matrices A,-, B;, and C,- are large sparse matrices. For all i stations 
except those on the airfoil surface (i.e., from i = ILE to i = ITE) each of blocks are 
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of size 3 J x 3 J, because there are three equations per node, for variables <f», x , and 
y. The variables are ordered so that 


~n+l 
X i,l 

7<T l + 1 

Vx,l 


-*"l 

-2/"i 


_ $71 
t 

x? +1 - X? 


_ (B? . 

T n +1 ™n 

?/m +1 - y?j 


(3.50) 


*”+ 1 _ , 
t 

~n+l n 

X i,J X i,J 

-, n +l — n n / 
Vi.J Vi,J / 


For i-stations on the airfoil surface, there are an additional two equations per 
node to compute the stagnation point location and grid movement along the airfoil 
surface, as described earlier. Hence, the blocks at these stations are formally of size 
(3 J + 2) X (3 J + 2) (although one of these equations is actually not necessary at the 
i = ILE station, since there is only one additional equation at the stagnation point). 
The arc length variables and Si : j are stored at the end of the vector given in 
Eq. (3.50) at these stations. 

The entries in the global matrix shown in Eq. (3.48) and its corresponding residual 
are computed in two ways. The terms due to the steady flow governing equation are 
assembled from the elemental stiffness matrix and force vector given in Eqs. (3.21) 
and (3.22), respectively. The stiffness matrix is evaluated at each computational cell, 
and its entries are assembled into the appropriate locations in the global matrix. The 
force vector is also evaluated at each computational cell, and its entries are assembled 
into the appropriate locations in the residual vector. 

The entries due to the grid generation equations, Eqs. (3.3) and (3.4), and the 
boundary conditions described in the previous sections are computed by linearizing 
the grid discretized equations and boundary conditions with respect to the value of 
the steady potential and grid location at the current iteration. For example, earlier 
in this chapter [Eq. (3.25)] we said that the boundary condition on the y-coordinate 
on the upper (j = J) surface of the computational domain on the periodic boundary 
may be written as 

Mf tJ = yi,j - yi, i - G = 0 (3.51) 

To perform Newton iteration on this equation, we rewrite it in the form 


dMjA 

dyi,j J 
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(3.52) 
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(3.53) 
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where the right-hand side of Eqs. (3.52) and (3.53) is the residual of Eq. (3.51) 
evaluated at the current iteration. The remaining boundary conditions and the grid 
equations may be linearized with respect to these variables in a similar fashion. Once 
the linearization has been performed, the entries need only to be assembled into their 
proper locations. The resulting linear system of equations may then be solved very 
efficiently using LU decomposition. 

After each Newton iteration, the grid locations and potential are updated to obtain 
a new estimate of the solution. This process is repeated until the error in the residual 
vector |M,Nj 

T is less than a specified tolerance. 

3.2 Unsteady Flow Solution Procedure 

Now that the steady solution has been obtained on a streamline grid, the unsteady 
flow may be calculated. For unsteady flows with vorticity (i.e., the forced response 
problem), the drift function and rotational velocity must be calculated at the before 
the assembly of the unsteady field equations and boundary conditions. For unsteady 
flows due to blade motion (the flutter problem), the unsteady grid motion must be 
calculated before the unsteady field equations and boundary conditions. 

3.2.1 Drift Function and Rotational Velocity Calculation 

Drift Function Evaluation 

The first step in the vortical analysis is to calculate the drift function, A, since the 
stream function, \H, is already known from the grid generation equations. Since the 
grid follows the flow streamlines, it is sufficient to define the stream function at each 
point along the inlet ( i = 1) boundary. From the discussion in Chapter 2 [Eq. (2.49)], 
the stream function may be expressed as 

df* j = R~ooU-co{y\,j yi,i) Q<x>{yi,j yi,i) (3.54) 

where R-oo and U-oo are the upstream steady density and x-velocity, the product of 
which is equal to the prescribed mass flux at the exit boundary of the computational 
domain, Q <x,. Note that the value of the stream function on the stagnation streamline 
on the lower boundary of the computational domain, iko, is assumed to be zero. 

For this analysis, the drift function will be set to zero at the upstream far-field 
boundary. Because the drift function measures the relative time it takes a fluid 
particle to move along a streamline, the initial value may be any arbitrary constant. 
Zero has been used here for convenience. Since the steady solution procedure resulted 
in a streamline grid, the calculation of the drift function at each grid location is 
relatively straightforward. Using the definition from Chapter 2 [Eq. (2.22)], the drift 
function may be calculated using the compound expression 


Ai'-ij + T (y«,i y*’— id) ]/ > 1 


(3.55) 
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The drift function is measured relative to the upstream far-field boundary, so the 
drift function at the i = 1 axial grid station is zero. Once the drift function has been 
calculated, the rotational velocity, v R , may be calculated. 


Rotational Velocity Calculation 

The calculation of the rotational velocity, v R , may be accomplished largely through 
Eqs. (2.37) and (2.54) described in Chapter 2, i.e., 


v *=K c ' + ai) VA + 


^2 + 


uy 




exp[j(KiA -f I < 2 ^)] 


(3.56) 


and 


7 3 \ i 3 R— oo oo G COS (C2-/^1 C1A.2) 

* = /d c,+ mttwT s ’ 
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2tt(4> - To) 

O P-—O0C COS 0 / QQ 


x exp[j(iGA + I < 2 ^) ] (3.57) 

Here the constants c x and c 2 , and the “wave numbers” K\ and K 2 may be calculated 
using the upstream divergence-free condition of v R [Eq. (2.33)] and the choice of the 
magnitude and phase of the vorticity. The stream function, '5, may be determined 
as part of the grid generation procedure, and the description of the drift function 
calculation was given earlier. All of the upstream flow conditions are prescribed. 
Hence, the only remaining tasks to determine v R are the calculation of the stagnation 
point constant a 0 and the gradients of the drift and stream functions. 

Recall from Eq. (2.28) that the constant a 0 is defined as 


0,0 — ~ 


a|v$|\ 

dn ) 


(3.58) 


Numerically, do may not be evaluated exactly at the stagnation point. Instead, we 
wish to calculate do at the midpoint between the stagnation point ( i = ILE ) and 
the previous adjacent grid point ( i = ILE — 1). Using a simple finite difference 
expression, ao may be written as 


do 


$ILE, 1 ~ $ILE- 1,1 1 ( SHE ~ SHE- 1 ) 

( SHE ~ SHE- 1) 2 


(3.59) 


where 

S[LE ~ SlLE- 1 = \J{*ILE,\ ~ ^ILE- 1,1 ) 2 + {VILE, 1 ~ VILE- 1,1 ) 2 (3.60) 

Numerical experiments have shown that v R is relatively insensitive to the details of 
the numerical calculation of ao, be., different numerical formulations of a 0 did not 
result in a measurable difference in the resulting unsteady potential. 

Finally, the gradients of the drift and stream functions must be evaluated. A 
general and straightforward approach to numerically calculating the gradient of a 
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function is to use a generalized coordinate transformation [54], Specifically, the gra- 
dient will be evaluated in the (E ,H) grid generation coordinate system. This may be 
accomplished by considering the derivatives with respect to x and y using the chain 


rule, i.e., 

d dZ d dH d 

dx dx dZ dx dH 

(3.61) 

and 

d dZ d dH d 
dy dy dZ dy dH 

(3.62) 


Ideally, however, we only wish to evaluate derivatives with respect to and H . This 
may be accomplished using the transformation matrix [J], which is defined as 


[J] = 


r dz 

dZ 1 

dx 

dy 

dH 

dH_ 

_ dx 

dy . 


(3.63) 


To obtain the derivatives we desire, the matrix [J] may be inverted, resulting in 


dZ dy/dH dH _ dy/d E 

dx ~ |[J]| dx ~ |[j]| 


(3.64) 


where 


dZ dx/dH dH dx/dZ 

Sy - |[j]l dy - |[J]| 


l[J]l = 


/ dx dy 
{ dZdH 


dx dy\ 1 
dHd ZJ 


(3.65) 

(3.66) 


Combining the above expressions, the derivatives in Eqs. 
rewritten as 


d_ 

dx 



dy_J_\ 
dZ dH) 


(3.61) and (3.62) may be 
(3.67) 


d irTll ( dx d dx d \ 
d^~' [ ^\dHdZ + dZdHJ 


(3.68) 


Note that now all derivatives are with respect to E and H. These derivatives may be 
computed using centered finite differences modified for variable spacing in E and H. 
The gradients of the drift and stream functions (or any other scalar) may be evaluated 
using this approach. Once the gradients have been calculated, the rotational velocity 
may be determined at each computational node. 


3.2.2 Unsteady Grid Generation 

For flutter problems, the motion of the computational grid must be computed before 
the calculation of the unsteady flow. The motion of the grid boundaries is easily 
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specified. The motion of the grid on the airfoil surface is required to match the 
vibrational mode shape of the blade motion. There is no grid motion in the far 
field to simplify the implementation of the far-field boundary conditions. Finally, the 
periodic and wake boundaries must smoothly connect these specified boundaries. 

Previous researchers have solved Laplace’s equation to determine the grid motion 
throughout the computational domain [30, 35, 24]. Laplace’s equation was used so 
that the grid motion varies smoothly throughout the computational domain. As a 
result, derivatives of the grid motion may be computed with minimal truncation error. 
Furthermore, in finite difference discretizations using centered differences, grid motion 
that does not vary smoothly may excite sawtooth modes in the unsteady solution. 

In the present method, sawtooth modes are not admitted into the finite element 
solution, so the grid motion need not be as smooth as for centered finite difference 
schemes. The truncation error of the derivatives must still be kept to a minimum, 
however. In this report, the calculation of the grid motion will be performed using 
a linear distribution of the grid motion over a few adjacent grid lines. A linear dis- 
tribution is quite simple to implement, results in a relatively smooth distribution of 
grid motion, and does not significantly increase the overall truncation error of the 
unsteady solution method. For most cases examined in the present work, approxi- 
mately five grid nodes in each direction were sufficient to obtain the correct unsteady 
solution. Mathematically, this may be expressed as 

L,j = Li Lj f A iR (3.69) 

where f A iR is the grid motion at the nearest grid point on the airfoil surface, and 
Li and Lj are functions that contain fractional values to distribute the grid motion 
smoothly. 

3.2.3 Finite Element Discretization 

The next step is to compute the stiffness matrix and force vector for the unsteady flow 
equations. The variational principle governing the unsteady flow was developed in 
Chapter 2 [Eq. (2.74)]. As for the steady flow, the variational principle is discretized 
using bilinear quadrilateral isoparametric elements. Since the unsteady flow problem 
is linear, however, the solution may be computed in one iteration. At each element 
n, the values of the unsteady potential at each of the corners are interpolated into 
the interior of the domain using the interpolation 

«{,-)) = |NJ.*„ (3.70) 

Similarly, the gradient operator may be written as 

V'0 = [N'U n (3.71) 

Substituting these expressions into the unsteady functional given in Eq. (2.74), results 
in 

IW = i /jf R ^ [0f[NTCvW* r [N'] n 0 n 
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+J'<» [NJIV'$ t [N'J^„) 

W^IlnJ I |NJ „</>„]} d(d n 

- [L R {V'$ T [j][N']„5„ + V' • fV'4 T [N'W„ - juV' ■ f |NJ„0„ 

- i i ? V'® 3 '[j]V'$V'$ T [N']„0 n + ^V'^[J]V'« LNj„?„ 

+ ^f • V'4>V'«V[N'W„ + Af . V'4>[Nj„?„ j dt dn 

+ [[ V • i?v R |NJ n </> n d£ drj (3.72) 

Setting the first variation of this functional to zero results in the discretized governing 
equations for the unsteady flow 

OT„„„ = j *{-{N'£[N'] n + i [[N'KV'<I»V'$ t [N']„ 

+j u ([NIJV'SLNJ. - |Nj;[V'4> T [N'] J1 ) + <o 2 |NJJ|NJ„]}0„ ii dr, 

+ JL r { yijfT ( lI] " 5l v ' ,i ' v ' $T ) l N 1" - ' V '*LNJ„ 

- (v'$ r [J] + V' ■ fV'<[> T - ^jjV'$ r [j)V'»V'# r ) [N']„ 

+;w [v' ■ f - J^V'$ T [j]V'$] LNj„} d( dr, 

+ if V ■ flv B |NJ„ d(dr,= 0 (3.73) 

J */£ n 

From this form of the first variation, it is clear that the local stiffness matrix, [k] n , 
may be written as 

[k]„ = / )f _ R {-(Nt[N']„ + i [[N'£VW$ T (N']„ 

+io.([N')JV'4LNJ„ - (NJ^V'$ t [N']„) + a, 2 |NJJ|NJ„]} d( dr, (3.74) 
and the corresponding unsteady force vector, e„, is 

e„ = / U |^f T ([I] - [N'] n - ■ V$|NJ„ 

- (V$ T [J] + V' • fV'$ T - ^V'$ r [J]VW$ r ) [N'] n 

+ju [v # . f - ^V'$ T [J]V'$] [NJnj + V • Rv r I NJ n ) d£ dr, (3.75) 

As in the steady flow calculation, the elemental stiffness matrix and force vector are 
computed for each cell in the computational domain and assembled into a global 
stiffness matrix and force vector. The next section describes the unsteady boundary 
conditions which also contribute terms to the global stiffness matrix and force vector. 
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3.2.4 Near-Field Boundary Conditions 

Periodic Boundary Condition 

In Chapter 2, it was noted that the periodic boundary condition would be enforced 
as part of the numerical solution procedure. This is not the only way to apply this 
condition, however. In Hall’s original linearized potential formulation [24], a Lagrange 
multiplier was used to enforce cascade periodicity on the upstream periodic boundary. 
Although this approach will result in the correct solution, it does not lend itself well 
to the sensitivity analysis procedure. Consequently, in this work a slightly different 
approach is used to enforce cascade periodicity. 

Examining the periodic boundary condition, Eq. (2. T9), shows that the unsteady 
potential on the upper and lower periodic surfaces of the computational domain are 
not independent. Hence, the value of the unsteady potential on one of the periodic 
boundaries need not be calculated explicitly. The number of equations at each axial 
station in the periodic region may then be reduced by one. This new vector of the 
discretized unsteady potential and the original vector are related by a simple linear 
relationship. If the original unsteady potential vector is denoted by </>, and the new 
reduced vector by </>, then 

<t> = [I] 4> (3.76) 

where 

(3.77) 


Note that the matrix [I] is rectangular, of size J x ( J — 1). 

The variational principle describing the flow may be written as 

n = -</> T [k\<f) — <p T e = min (3.78) 

where [k] is the global stiffness matrix and e is the force vector. Substituting the 
solution vector transformation equation, Eq. (3.76), results in 

n = |$ r f7r{k][/j5 - J T (7] r e (3.79) 

where the overbars again denote the complex conjugate. Setting the first variation of 
this expression to zero gives 

SI I = [7] T [k] [T\<j> - \j] T e = 0 (3.80) 

This is the equation that must be computed. In practice, Eq. (3.80) means that the 
periodicity condition is enforced by pre- and/or post-multiplying the stiffness matrices 
and force vectors by [I] r and [I], respectively, once the block stiffness matrices and 
force vectors have been computed. Hence, the block stiffness matrices [a,], [b,], and 
[cj] and block force vector e t at each station i are replaced according to the following 
table 


ffl = 
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i — 

1 

2 ILE-1 

ILE 

ILE+1 

ILE+2 -> . . . 

M -»• 

— 


[I] T [a,][i] 


w 

[bi] — > 



[i] T N[i] 

[b,] 

[bj 

M -> 

[i] T [c t ][i] 

[i] T [c,][i] 

ma] 

N 

N 

e,- -*• 

[I] T e, 

ft 3 * 

ft 3 * 

e. 

e, 


So, for example, from axial grid stations i = 2 to i = ILE — 1, the block equation to 
be solved is of the form 

1 + [b,\4> t + [c i\<f> l+1 = e, (3.81) 

where [a,], [b,], [c,], and e, are the block stiffness matrices and force vector, respec- 
tively, that have been modified according to the table given above. Although this 
may appear to be somewhat complicated, the sparseness of the transformation ma- 
trix actually results in very few changes to the computed stiffness matrices and force 
vector, and eliminates the need to use Lagrange multipliers. 


Airfoil Surface 


The boundary condition on the airfoil surface is that there is an upwash on the surface 
of the airfoil due to blade motion and if there is a nonzero rotational velocity on the 
blade surface. Recall from Chapter 2, that the boundary integral terms in the first 
variation of the unsteady variational principle may be written as 


linear = J J^{- ■ •} H ^7? 

+ j> R (joj f • n - [J]V$ • n + V K ■ n - ^ 


6<f>ds = 0 


(3.82) 


The first two terms in the boundary integral arise as the natural boundary condition 
of the variational principle. The third term, however, was added separately to account 
for the upwash due to a nonzero rotational velocity on the airfoil surface. The integral 
of this third term is discretized using one-dimensional finite elements, and the result 
is added to the inhomogeneous term e. On the (j = 1) surface, the equation to be 
solved is 

R 

[k] = e,,i -| — ^ i — 2/»-ii0 T v y ,t-i,i( a: «,i — x *'-i,i)J 


Ri, 1 


— 2/»-l,l) T ( x i,l x i- 1,1 ) 


Ri. 




3/i,l) + w y ,*,l( x «+l,l X »,l) 
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+ 


-^+1,1 


- u £i+l,l(!fc+M - Vi, 0 + V y,i+l,l( x i+hl ~ *i.l)] 


(3.83) 


where the x and y subscripts refer to the x- and y-components of the rotational 
velocity v R , respectively. On the (j = J) surface, 


Rt 1 

= e i,J g ~ Vi- 1,1 ) d" V y,i-l,l( X i,l ~ ^i-l.l)] 


[~ v x,i,i (y*M - yi-i.i) + ^,1(^,1 - 

D 

[-«&•, ito+u - j/i,i) + i&itei+i.i - *< t i)] 
~ V x,i+l,l(y*+lA ~~ y*.l) d" U v,i+l,l( X *+l,l — x i, l) 


3 

6 


(3.84) 


For airfoils containing a stagnation point, of course, it was shown in Chapter 2 that 
v ft is zero on airfoil and wake surfaces using the Atassi modification. So this term 
need only be computed for airfoils that do not contain a stagnation point and the 
Atassi formulation is not used. 


Downstream Wake Boundary Condition 

As was shown in Chapter 2, there are two separate parts to the wake boundary 
condition that need to be considered. The first is the auxiliary pressure continuity 
condition, Eq. (2.86). In discretized form, the pressure continuity condition for node 
i may be written as 


jwRi-i/2 ( Si 


Si— 1 ) 


(<fc- 1,1 d- <f>i, 1 ) 
2 


(j-i,J d- ^;,j) e 
2 


d" Ri-XllVi-X^ (Si — 5j-l) 


(0»,1 ~ l,l) 

(s,' ) 


(•Si - 3,-0 


= 0 


where 




Si - 1 = y/(xi - x^ 0 2 + (yi ~ Vi- 0 2 


(3.85) 

(3.86) 


Equation (3.85) provides closure for the wake displacement variable r;. 

The other condition is that mass must be conserved on the wake surface. Specif- 
ically, as shown by Eq. (2.88), there is an additional upwash that must be added 
due to the displacement of the wake. This additional upwash is added on both wake 
surfaces in the computational domain. On the ( j = 1) surface, the equation that 
must be solved is 

[k]i,i fa - e tll 


--Ri+i /2 Vi+i /2 (r i+ i - ri) - -jwRi+y 2 (s i+ i - Si) (r t+a + 2r<) 
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- e (r,_! - 2 n + r i+1 ) = 0 (3.87) 

and 

•Si+i - 5,- = y(x,+i - Xi ) 2 + (y,+i - yi) 2 (3.88) 

The last term on the left-hand side of Eq. (3.87) represents second-difference smooth- 
ing applied to the wake displacement to prevent sawtooth modes. For most cases, 
the coefficient, e, is equal to 0.1. 

Similarly, on the (j = J) surface, 

- e t)1 

+ ^Ri + !/2 K+ 1/2 (r<+i - r t ) + ^juRi+i /2 (s i+ i - Si) (r , +1 + 2 r<) 

+ -Ri-i /2 Vi- 1/2 (r; - r»-i) + gj'w-R,-- 1/2 (s,- - s t _! ) (2r,- + r^) e- 7 * 

+ £ (r,_! - 2r t - + r t+1 ) e J<r = 0 (3.89) 

With the addition of the wake displacement variable, r,-, and the boundary con- 
ditions given above, note that there are now J + 1 equations for J + 1 variables at 
each i-station in the wake. It is useful here to define a new vector to represent the 
combination of the potential at each node and the wake displacement, i.e., 

+* - (-£-) ( 3J °) 

Hence, including the auxiliary boundary condition for the wake displacement, the 
block stiffness matrices [a*], [b,], and [c,] and block force vector e, are replaced ac- 
cording to the following table 


i = 

. . . -> ITE-2 

ITE-1 

ITE 

ITE+1 -*■ IMAX-1 

IMAX 

[a.-] -*■ 

W 

[a»] 

M 

[a,] 

[^] 

M -*■ 

N 

N 

[b,3 

[bi] 

[b,] 

M- 

N 

[c»] 

N 

fe] 

— 

e, -> 

e. 

e; 

o 

ei 

e t 

e. 


So, for example, from axial grid stations i = IT E + 1 to i = I MAX — 1, the block 
equation to be solved is of the form 

[a,j<k-i + [bt]0< + [ctM’+i = e t - (3.91) 

where [a,], [b t ], [c,], and e,- are the block stiffness matrices and force vector, respec- 
tively, that have been modified according to the table given above. 

Finally, the computation of the rotational velocity upwash is calculated on the 
wake surface the same way as on the airfoil surface, using Eqs. (3.83)- (3.84). 
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3.2.5 Far-Field Boundary Conditions 

In Chapter 2, the basic analytical description of the unsteady flow in the far field was 
developed. In this section, the numerical solution of these equations is considered, 
including the presence of unsteady vorticity. In both the upstream and downstream 
far-held regions, the steady how is assumed to be uniform. If, in addition, the com- 
putational grid is regularly spaced in the axial direction, the unsteady how at each 
axial grid station i satishes the equation 


[»,]</»,_! + [b + [c t ]<£ 1+1 = e, (3.92) 

where the inhomogeneous term e, results from either non- divergence-free unsteady 
vorticity in the how, or a specihed acoustic gust. This inhomogeneous equation may 
be solved by examining the combination of the homogeneous part of the potential, 
4>^ , and the particular part of the potential, </>f. The goal of this analysis is to com- 
pute these two potentials so that the combination of the two results in the unsteady 
potential that makes the far-held boundaries nonrehective, i.e., 

& = 4? + 4f (3.93) 

Substituting this expression into Eq. (3.92) results in 

M (0?-, + 0f-O + M (0? + 0f) + M (0", + 0f +1 ) = e. (3.94) 
which may be written as two separate equations: a homogeneous equation 

WA + [bi)0f + [c i ]0't, = 0 (3.95) 

and an inhomogeneous equation 

+ [b i]4>f + [c = e t (3.96) 


Homogeneous Solution 

The solution of the homogeneous problem may be computed by considering how 
pressure and vorticity waves propagate from one axial grid station to the next. For 
example, if represents the homogeneous potential at grid station i, how can 
be represented? The desired relationship between the two is of the form 

0."-! = [T] 0" (3.97) 

where [T] is a matrix that acts to hlter out pressure waves approaching from outside 
of the computational domain. This matrix may be constructed in three interme- 
diate steps. First, the as yet unknown potential vector at the far held is Fourier 
transformed, so that 

1 [G 

di,m = 77 / <j>{xi,y)exv{-jPrny) dy 
Lr Jo 


(3.98) 
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where d, im are the Fourier coefficients at axial grid station i. In the present method, 
the coefficients are computed using the trapezoidal rule. Typically, the coefficient 
vector d; is much smaller in length than the number of computational nodes in the 
circumferential direction, since truncation error increases greatly for higher harmonics 
of the principal circumferential wave number. For most cases considered here, five 
Fourier modes (—2 < m < 2) are sufficient to describe the unsteady flow. In matrix 
notation, Eq. (3.98) may be written in the form 

d, = [Fxj^f (3.99) 

where [F a ] is a (generally rectangular) matrix that represents the discretized form of 
the trapezoidal rule. 

The next step in the analysis is to filter and propagate these Fourier modes to the 
next grid station. In the far field, since the steady flow is assumed to be uniform, the 
wave numbers of the wave propagation are known analytically. From the analysis in 
Chapter 2, the harmonic behavior of the potential may be written as 

OO 

<f>(x,y) = J2 d m exp[ja m (x - x ff ) + jf3 m y] (3.100) 

771 = — OO 


where is the x-coordinate of the far-field boundary being examined. Therefore, if 
the spacing between nodes on adjacent axial grid stations is known, then the propa- 
gation of the Fourier coefficients may be expressed as 

di- i,m = d i>m exp [~ja m (xi - x,_i) - j/? m (y t - y;_i)] (3.101) 

where 

U( to + (3 m V)± Cy/0l(U* + V 2 - C 2 ) + + 

arn ~ (J2 _ JJ2 

and 

A. = (3.103) 

In matrix notation, the filter and propagation relationship may be written as 

d,-i = [F 2 ]d, (3.104) 

where [F 2 ] is square and diagonal, and the diagonal entries are the above exponential 
expression for outgoing modes; the entries corresponding to incoming modes are zero 
'so that their contribution is not included in the final solution. 

The third step in the solution process is to perform an inverse Fourier transform 
on the Fourier coefficients d,_i to recover the the new values of the potential at each 
node. In matrix form, this may be written as 

<j>”_ x = [F 3 ] d,_ a (3.105) 

where the (also generally rectangular) matrix [F 3 ] is the discrete matrix form of the 
inverse Fourier transform. 
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Taken together, these three steps result in 

= IF,) [F,] [F.l <A" = [T] (3.106) 

so the transition matrix, [T] , is now known. Now that the behavior of the continuous 
part of the potential has been described, the particular part due to incoming gusts 
may be calculated. 


Particular Solution 


The particular solution may be calculated by considering the propagation character- 
istics of the incoming acoustic or vortical gust. Because the waves convect with the 
steady flow and the steady flow is uniform in the far field, the convection relation 
may be expressed as 

<t>i+n = 4>i (3.107) 

where 

= exp [jna p ( x i+1 - X;) + jn(3 p (y i+l - y { )] (3.108) 

For a pressure gust, the wave numbers a p and /? p have the same form as those for 
ct m and /3 m as shown above in Eqs. (3.102) and (3.103). For a vortical gust, the wave 
numbers a v and {3 V are defined as 


a 


v 


+ P P V) 

u 


(3.109) 


and 

ft = § (3.110) 

with U and V being the x- and .(/-components of the steady velocity at the far field. 

Consequently, the equation that describes the particular solution in the far field, 
Eq. (3.96), may be written as 


1 

z 


+ M*f + Z [c,)0f = e. 


(3,111) 


So the particular solution may be computed by solving the above linear system, 


i.e.. 



~ z [ a .] + [bt] + z [c,] 




e t 


(3.112) 


Equation (3.112) gives a relatively simple result that is best applied when a vortical 
gust is present. For an acoustic gust, there is a simpler approach. The unsteady 
pressure is linearly related to the potential by the expression 


p{x,y) = -R\ju + V$ • V] <f)(x«, y s ) exp \ja p (x - x ff ) +j/3 p {y - y ff )] (3.113) 

where (xff ,yg) refers to a coordinate at the far-field boundary. Hence, the particular 
solution may be specified once the magnitude and phase of the gust has been chosen, 
i.e., 

4 > p (x{{, y s ) = r . exp ti a p( x ~ Xff ) + 7^(2/ - ys)] (3.114) 

-R\ju + j.a p U + JP P V\ 
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Now that the general procedure for calculating the homogeneous and particular 
parts of the unsteady potential have been described, the actual implementation for 
each boundary may be examined. 

Upstream Far-Field Boundary 

At the upstream far-held boundary, the homogeneous part of the potential satisfies 
Eq. (3.95) with i = 1. 

[a-i]$o + [bi]^j + [ci]0 2 = 0 (3.115) 

where the block matrices and potential vectors have been modified in the periodic 
boundary conditions using Eq. (3.76). The transition matrix from Eq. (3.106), may 
be modified in a similar fashion, i.e. , 

[ t ] - mm = [ t ] 

Hence, Eq. (3.115) may be written as 

([a,][T] + [b 1 ])5f+[c 1 ]$f = 0 

and in light of Eq. (3.93), some algebra results in 

[b*]<£i + [cj]^ = ([bj] + 2 [ci]) = el 

where 

[6;] = ([ai][T) + [b,]) 

Hence, the far-held boundary conditions modify the block matrix 
force vector €»i. 

Downstream Far-Field Boundary 

The downstream far-held boundary conditions are assembled in an analogous fashion 
at axial grid station I, with the appropriate change in the direction of propagation of 
the pressure modes in the homogeneous solution. The only major difference is that 
the wake must be accounted for. Recall from Chapter 2 that there is a continuous 
and discontinuous part to the unsteady potential in the downstream far held 

<t> = <j* c + <t> d (3.120) 

Both the continuous and discontinuous parts of the potential must propagate cor- 
rectly. The transition matrix, [T], derived earlier, only applies to the continuous part 
of the potential. The discontinuous part propagates out of the domain as a vorticity 
wave. Hence, a new transition matrix may be dehned as 

[T*] = [T][I — D] + z[D\ (3.121) 

where [D] is the matrix operator that computes the discontinuous part of the poten- 
tial, and [I] is the identity matrix. The J + 1 equation is the propagated form of the 


(3.116) 

(3.117) 

(3.118) 

(3.119) 
[bi] and the block 
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mass conservation boundary condition. Since there is no grid motion in the far field, 
the conservation of mass boundary condition is 


d(j) d$ dr 

— + — — 

In discrete form, Eq. (3.122) may be written as 

9 ± ju> , , , d<t> (r /+1 - r/) 

■x~d>/ “ ~ 7 C ( r /+i + ri) - -r-7 r = 0 

dn 2 ds (si — Si-i) 


(3.122) 


(3.123) 


Solving for r/ +1 gives 

9 /dn -ju/ 2 + (d$/ds)/(sj - Si- 1) _ 

;w/2 + («*/&)/(*- *-i) ' >u/2 + (0*/3»)/(*-»i-,) r ‘ r,+1 V' lM > 


or 


L^J <£/ + £r/=r /+1 (3.125) 

So the new transition matrix, with the additional row and column for the wake 
displacement, has the structure 


/ 


\ 


4>i+i 

\ r/ + i / 


[T] 


0 




(3.126) 


V r/ / 


In the wake boundary conditions, we introduced a new vector that contains the po- 
tential and wake displacement, so that 


4>i = 




Ti 


Hence, Eq. (3.126) may be written as 

4>»i = [TJ0, 


(3.127) 


(3.128) 


Using this new transition matrix, the matrix equation at the downstream far field 
may be written as 


where in this case 


[&,]&_, + [VM>i = Q n + [£;]) if = 

|b;] = ([b,] + [c,][T]) 


(3.129) 

(3.130) 
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3.2.6 Assembly and Solution 

Like the steady flow equations, the unsteady flow equations may be assembled into 
a block-tridiagonal form. Since there is only one variable per node to be computed, 
however, the size of the blocks of the matrix is significantly less than the steady 
flow case. Schematically, the matrix structure may be written using the notation 
developed in this chapter as 

b* Ci 

5 2 b 2 c 2 

a,' b, C{ 

a/_! b/_i c/_j 

a/ 6 J 

or, more generally, 

[W]{0} = {e} (3.132) 

In the periodic region, the blocks of the matrix are of size (J — 1) x (J — 1), since 
the periodic boundary condition eliminates one equation from each block equation. 
For i-stations on the airfoil surface [signified by equation i in Eq. (3.131)], the block 
matrices are of size J x J, because the boundary conditions do not change the number 
of equations in this region. In the wake region, the block matrices are of size ( J + 
1) x (J + 1) due to the additional wake displacement variable and its corresponding 
equation. Recall that the asterisk indicates that the block matrix or right-hand side 
vector are modified due to the application of the far-field boundary conditions. 

The entries in the block matrices and right-hand side vectors are taken directly 
from the stiffness matrix and force vector computed at each computational cell with 
additional boundary conditions added in the appropriate locations as described in 
this chapter. Once assembled, this linear system may be solved very efficiently using 
LU decomposition to determine the unsteady potential </> at each computational node. 


01 


e* 

c i 

02 


e 2 

■ •• 

► = < 

e, 

0/-i 


ej_i 

4>r 


- - 


(3.131) 



Chapter 4 

Sensitivity Analysis 


Chapters 2 and 3 described a general method to calculate the steady and small dis- 
turbance unsteady flow over a two-dimensional cascade of airfoils due to blade motion 
and incoming gusts. The next task (and the primary goal of this report) is to deter- 
mine the effect small changes in the airfoil shape or cascade geometry have on the 
steady and unsteady flow fields. In Section 4.1, the general procedure for this sensi- 
tivity analysis will be described. The sensitivity analysis makes use of the nominal 
steady and unsteady flow LU decompositions so that no additional matrices need 
to be factored, resulting in a very computationally efficient procedure. Section 4.2 
contains a description of the steady flow sensitivity analysis procedure, including a 
discussion of how the perturbation of the airfoil shape is prescribed. In a similar 
fashion, the sensitivity of the unsteady flow is considered in Section 4.3. 

4.1 General Procedure 

The quantities of interest to be calculated in this report may be separated into four 
groups: first is the nominal steady potential and grid locations, $ and x; second, 
the perturbation of the steady potential and grid due to small changes in the airfoil 
geometry, and x'; third, the nominal unsteady potential, <f>] and fourth, the per- 
turbation of the unsteady potential due to small changes in the airfoil geometry, <fi' . 
Chapters 2 and 3 described how the nominal steady flow (i.e., $ and x) and nominal 
unsteady flow (i.e., 4>) are calculated. This section describes a general procedure for 
calculating the remaining quantities, which are referred to as perturbations of their 
associated nominal quantities. Throughout this chapter, a prime (') will be used to 
denote these perturbations. 

It should be noted that the perturbation of the steady potential, and the 
nominal unsteady potential, <fi, are not of the same order. Figure 4.1 illustrates these 
two separate perturbations. In Chapters 2 and 3, the perturbation series we examined 
were expanded using a small parameter on the order of the blade displacement, f. 
Here we expand the equations in perturbation series using a small parameter on the 
order of the change in geometry, which is denoted here by e. Hence, the perturbation 
of the unsteady potential, cf)' is not formally second order. Instead, it is first order in 


76 



4.1. GENERAL PROCEDURE 


77 


Nonlinear Flow Equations 


Expand in small parameter (blade displacement) to get: 


Nonlinear Steady Flow 


Linearized Unsteady Flow 

Discretize to obtain: 


Discrete Steady Eqns. 

0(1) 


Discrete Unsteady Eqns. 

o( f) 


Expand in small parameter (change in geometry) to get: 


Steady Sensitivity 


Unsteady Sensitivity 

0(e) 


0(e- f) 


Figure 4.1: Overview of quantities in sensitivity analysis. 


both e and f . 

In Chapter 3, it was shown that the steady flow on a streamline grid may be com- 
puted by solving a set of nonlinear equations, Eqs. (3.8) and (3.9), for the unknown 
grid position and steady potential. Taken together, these equations may be written 
as 


N(*,x;Z) 

M(*,x;Z> 


(4.1) 


Because this set of equations are nonlinear, we chose to use Newton iteration to reduce 
the set of nonlinear equations to a series of linear equations to solve. The Newton 
iteration procedure [Eq. (3.47)] may be written in the form 


1 $ 

. n+1 

f $ \ 

n 

r <9N 

9$ 

<9N 1 
dx 

V x . 

) = 

l X ) 


§ CC5 
1 

dM 
dx -1 



(4.2) 


When the residual of the iteration procedure has been reduced to a specified tolerance, 
the nominal steady potential and grid locations are obtained. 

This section will consider the effect a small change in geometry has on these 
solutions. Consider a small perturbation Z' in the vector describing the airfoil shape 
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and cascade geometry. This will result in a perturbation of the steady potential and 
the grid location. The perturbed solution will satisfy the equations 

N($ + $ / ,x + x , ;Z + Z') \ 

M($ + $',x + x';Z + Z') J 


where the primed quantities are small perturbation quantities arising from the pre- 
scribed perturbation in the geometry Z' . Expanding Eq. (4.3) in a Taylor series about 
the nominal solution gives 


r <9N 

<9N I 


dx 

dM 

dM 

L a$ 

5x J 


/ £N Z , \ 

~dZ Zl 


dM 7 , 
~ dZ Z 


(4.4) 


where we have grouped the terms so that left-hand side of Eq. (4.4) contains terms 
involving the perturbations of the dependent variables # and x, and the inhomo- 
geneous right-hand side contains terms involving the prescribed perturbation of the 
independent variable Z. To solve for the unknown perturbations $' and x', Eq. (4.4) 
is rearranged slightly to obtain 


/ \ 

U'/ 


r dN dN 1 

-l 

( 5N 7' \ 

a$ dx 


dz L 

dM dM 


dM„, 

ax J 


dz z / 


(4,5) 


Note the similarity of Eq. (4.5) to Eq. (4.2). The same matrix must be “inverted” to 
obtain the perturbed steady solution that was used in the last iteration of the Newton 
solver. In other words, the same system of equations as in the nominal calculation 
must be solved, only with a different right-hand side due to the perturbation Z'. 
Therefore, if Eq. (4.2) has been computed using LU decomposition, and the factored 
matrix has been saved, then the sensitivity analysis may be obtained for virtually no 
additional computational work. This is because the perturbations and x' may be 
computed through a simple back-substitution, which requires significantly less com- 
putational time than a complete solution. Furthermore, the relative computational 
savings increases with the number of design variables examined. The perturbation 
of each design variable will result in a new right-hand side of Eq. (4.5), each only 
requiring a back-substitution, so many design variables may be examined before the 
cost of the sensitivity analysis approaches the cost of the single nominal calculation. 

Consequently, the main task of the steady flow sensitivity analysis is to determine 
the prescribed perturbation Z' and its associated matrices dN jdZ and dM/dZ. The 
details of the computation of these quantities will be examined in Section 4.2. 

Having computed the perturbation of the steady potential and grid due to a 
change in geometry in addition to the nominal solution, we now consider the small 
disturbance unsteady flow. From Chapter 3, we know that the discretized unsteady 
flow may be represented by a linear system of equations of the form 


[W(*,x,w,<r)]{^} = {e($,x,v H ,f ,u,<r)} 


(4,6) 
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After solving Eq. (4.6), the perturbation of the unsteady potential due to a change 
in the geometry may be calculated. The governing equations of the perturbed un- 
steady flowfield due to small changes in geometry, frequency, rotational velocity, in- 
terblade phase angle, and mode shape will be of the form 


[W($ + x + x', w + u/, a + a')] {cf> + 0' } 

= e(3> + $',x + x', v R + v R ',f + i',u + u>',a + a 1 ) (4.7) 

Expanding the discretized form of Eq. (4.7) in a perturbation series and collecting 
terms of first order gives the desired equation for the unknown discrete approximation 
to the perturbation of the unsteady potential, <//, i.e., 


[W]0' = 


de 


$' + 


de 

!ht 


x' + 


de 

dv R 


v R ' + 


'de 

df 


f' + 


de 

du 


u' + 


de 

da 


dW 


<9$ 


+ 


dW 


dx 


x' + 


aw 


du 


J + 


aw 


da 


4 > 


(4.8) 


or, more succinctly 

[W]<£' = e' - [W]'(f> (4.9) 

Note that the computation of x' and has already been described. Section 4.3 
contains a discussion of the computation of u>' and f' for the flutter problem, and \ R ' 
for the forced response problem. In addition, Section 4.3 will contain a discussion of 
the assembly of the right-hand side of Eq. (4.8). 

Note that the sensitivity analysis procedure outlined in this section is not limited 
to the flow model and discretization scheme chosen in the present method. Equa- 
tions (4.1)-(4.9) may be applied to any linearized flow model, such as the linearized 
Euler or Navier-Stokes equations. Finite difference and finite volume discretizations 
also may be used. The only requirement to obtain the computational efficiency of the 
procedure described in this chapter is that the nominal flow equations must be solved 
using LU decomposition or some other algorithm designed to solve a linear system of 
equations for multiple right-hand sides. 


4.2 Steady Flow Sensitivity Analysis 

4.2.1 Prescribing the Perturbation 

The first step in the steady flow sensitivity analysis is to prescribe the perturbation 
in the vector containing the airfoil shape and cascade geometry, Z'. In general, 
the basis functions chosen to prescribe the geometry perturbation may be chosen 
arbitrarily. There is no limit on the number of functions that may be chosen, nor is 
there a restriction (such as orthogonality) on the choices. To provide physical insight, 
however, these functions should be chosen so that they are meaningful to the designer. 
In this work, the basis functions may be classified into two groups, cascade geometry 
and airfoil shape basis functions. 
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Recall from Chapter 3 that cascade geometry parameters refer to variables that 
define the cascade behavior independently of the actual airfoil shape. For example, 
the stagger angle, 0, and the blade-to-blade gap, G, are both variables that are 
independent of the airfoil shape. Airfoil shape parameters refer to variables that 
define the airfoil shape, for example, the magnitude of the camber and thickness. 

The details of the generation of the actual perturbation of the airfoil shape is 
highly dependent on the type of airfoil definition used. Although the details will not 
be discussed here, Appendix B contains a detailed description of how the perturbation 
of shapes from a common airfoil series (the NACA modified four digit series) may 
be generated. In the current discussion, however, it is sufficient to assume that 
perturbation of the airfoil shape or cascade geometry is prescribed and is contained 
in the vector Z'. 

Once the airfoil shape perturbation has been defined, the next task in the steady 
sensitivity analysis is to define the entries in the 3M/5Z and dN/dZ matrices shown 
in Eq. (4.4). These matrices are extremely sparse because nearly all of the inhomo- 
geneous terms in Eq. (4.4) arise from the boundary conditions. In practice, these 
matrices need not actually be formed; it is sufficiently straightforward to compute 
the matrix-vector product [the right-hand side of Eq. (4.4)] directly. 

There are three main sources of inhomogeneous terms in Eq. (4.4). First, the 
near-field boundary conditions, contribute terms because these conditions require the 
grid to remain on the airfoil surface. Any change in the airfoil shape, then, will 
influence the location of the grid points. Second, there is a contribution from the far- 
field boundary conditions primarily due to changes in the blade-to-blade gap. Finally, 
there are inhomogeneous terms due to the the grid equations, Eq. (3.8), because of 
the requirement that the grid follow the flow streamlines. 


4.2.2 Near-Field Boundary Conditions 


Periodic Boundary Condition 


Consider the steady periodic boundary conditions, Eqs. (3 .23)— (3.25) . Taking the 
derivative of the boundary condition equations with respect to the shape variables 
described earlier and multiplying by the perturbation of the vector containing the 
airfoil shape results in terms that contribute to the right-hand side of Eq. (4.4). 
Assuming that the fractional arc length array F % is not a design variable, the boundary 
conditions on the j = 1 periodic boundary of the computational domain do not result 
in inhomogeneous terms in Eq. (4.4), i.e, 


(4,10) 


am 

~3Z~ Z “ ~dZ Z ~ 

On the j = J surface, although the boundary condition governing the ^-location is 
only dependent on the grid, the y-location is also dependent on the blade-to-blade 
gap, so there is an inhomogeneous term associated with the boundary condition on 
the y-grid equation, i.e., 

8MIj 


dZ 


Z' = 0 


(4.11) 
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and 


dMl 

dZ 


Z' = 


-a 


(4.12) 


where G' is the perturbation in the blade-to-blade gap. Hence, the periodic boundary 
conditions are not directly dependent on the airfoil shape. 


Airfoil Surface Boundary Conditions 


The next set of boundary conditions to consider are those that constrain the grid 
points to lie on the airfoil surface, Eqs. (3.27)-(3.30). On the j = 1 surface, there are 
inhomogeneous terms that result from the spline evaluation, i.e., 


and 


dM ? i 
dZ 

mi 

dZ 


Z' = 


z' = 


(4.13) 


(4.14) 


where X' s and Y$ are the perturbations in the cubic spline evaluation of the airfoil 
surface coordinates corresponding to the arc length s,^. Similarly, on the j = J 
surface, 


dM f 


and 


dZ 

m,j 

dZ 


V = -X'sM 


v = -Yfa itJ ) - a 


(4.15) 

(4.16) 


where the latter equation also includes perturbations of the blade-to-blade gap, G' . 

The next two boundary conditions are those that require the grid to move propor- 
tionally to the movement of the stagnation point, Eqs. (3.33)-(3.34). These conditions 
are applied at every airfoil surface grid location except the stagnation point. On the 
j — 1 surface, 


dM\ 


t,J+i 


dZ 


Z' = 


while on the j — J surface, 


Wm r 
— W z - 


- c,_i a , 
1 - Ei-i 1 

(4.17) 

Fi — Fi- 1 c/ 
1 — Fi-i N 

(4,18) 


Recall that S\ and Sn are the arc length of the first and last points, respectively, of 
the set of points defining the surface of the airfoil, as shown in Fig. 3.4. 

The final surface boundary condition to consider is the condition at the stagnation 
point, Eq. (3.35). This boundary condition requires that the grid node at the leading 
edge of the airfoil be located at the stagnation point. The right-hand side term is 
found to be 


dM 


ILE.J+l 

dZ 


Z' 


1 


SILE+ 1,1 - s ILE+l,J 


S ILE S ILE+\,J vr n , 
\_ s ILE+\,\ S ILE 


(4.19) 
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The term shown here arises from the steady periodicity condition. The potential 
on the upper surface must be evaluated on the reference airfoil. From the steady 
upstream far-held boundary conditions it can be shown that the steady potential on 
the upper periodic surface and the steady potential on the lower periodic surface are 
related through the parameters G and V-oo. Specifically, 

$i,i = - G y_oo (4.20) 

so clearly there is an inhomogeneous term arising from perturbations in the blade-to- 
blade gap. 


Downstream Wake Boundary Conditions 

The final set of near-held boundary conditions to consider are those on the down- 
stream wake boundary, Eqs. (3.36)— (3.38). Recall from Chapter 3 that the steady 
wake boundary conditions are nearly identical to the boundary conditions applied on 
the periodic surfaces. Hence, like the periodic boundary, there are no inhomogeneous 
terms in the perturbation of the boundary conditions on the j = 1 surface, i.e., 


dM? x dM? x 

i# z =i ^ z = 0 


(4.21) 


On the j = J surface, the inhomogeneous terms account for perturbations in the 
blade-to-blade gap as well as the possibility of a blunt trailing edge, so that 


dMf t j 

dZ 


ZJ 


= 0 


(4.22) 


dMl 

dZ 


Z' = 


- [G' - (y; - j- )] 


(4.23) 


4.2.3 Far-Field Boundary Conditions 


Upstream Far-Field Boundary Conditions 


On the inflow boundary, only the equations governing the computational grid location, 
Eqs. (3.41)-(3.42), may result in inhomogeneous terms. Since the x-coordinate of the 
upstream boundary is hxed in this analysis, however, it is not considered to be a design 
variable, and therefore only the boundary condition on the y-coordinate actually has 
an inhomogeneous part to consider. Linearization of the boundary conditions results 
in the expressions 

"St z- = 0 

dZ 

and dMf . 

Z' = —FjG (4.25) 


(4.24) 


where it is assumed that the fractional arc length array in the far field, Fj, is not a 
design variable, i.e., it is not subject to change in the design process. 
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Downstream Far-Field Boundary Conditions 

The boundary conditions at the exit boundary are more complicated, mainly due to 
the possibility of a blunt trailing edge. It is assumed here that any perturbation in 
blade-to-blade gap G' is constant throughout the computational domain. Hence, if 
the airfoil being analyzed has a blunt trailing edge, the space in the grid representing 
the viscous wake from the trailing edge of the airfoil to the far-field boundary (see 
Fig. 3.5) is constant regardless of the perturbation of the gap. The size of the space 
representing the wake is only affected by perturbations in the first and last points of 
the airfoil definition. Consequently, the inhomogeneous terms arising from the flow 
boundary conditions, Eqs. (3.43)-(3.44), may be written as 

7 , . dNjj+i 1 G'(yi t j + 1 - yjj) 

dZ 8Z 2 ^°° G-{Y 1 -Y n ) 

+ (426) 

The terms arising from the grid boundary conditions, Eqs. (3.45)-(3.46), are similar 
to those for the inflow boundary, modified for the possibility of airfoil shapes with a 
blunt trailing edge, i.e., 


' = -Fj[G'-{Y{-Y^)} (4.28) 




(4.27) 


4.2.4 Sensitivity of the Stream Function 

The final step in the steady sensitivity analysis is to consider the inhomogeneous 
terms due to a perturbation in the stream function, \k' . Inhomogeneous terms arise 
because in this analysis the stream function at any node j is linearly related to the 
blade-to-blade gap, i.e., 

4b = QooFjG (4.29) 

where Q <*, is the prescribed mass flux at the exit of the computational domain and 
Fj is the fractional arc length array value at node j. Consequently, a perturbation in 
the gap results in a perturbation of the stream function, so that 


% = QooFjG' 


(4.30) 


Hence, any perturbation in the gap results in inhomogeneous terms in Eq. (4.4) arising 
from the field grid generation equations, Eqs. (3.3)-(3.4). These inhomogeneous terms 
may be written as 


^x 

~sz z = 


— / d 2 x — d 2 x _ d 2 x 

dEd-% ~ 2/j dEd¥ + 7 (^ 2 )' 


/ - 2 \ >dx —2 dx T 2 ^ dx 


(4,31) 
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and 


dMN 


,&y 

dl 


hi v' — - -> “ y ns' d 2 y r,-Q ” ~ . — 

Z a q -2 ^ dZd^ ^ 




dZ dE 2 ^ dEdV ^ dZdV 1 (M 2 )' 

+ (f)'Q^ + rQ'^ + fQ*L 


(4.32) 


where the primes in the denominator of some of the terms refer to the perturbation 
of the finite difference representation of the d/d $ operator. Note that the terms in 
Eqs. (4.31) and (4.32) only appear if there is a perturbation in the gap. 


4.3 Unsteady Flow Sensitivity Analysis 


4.3.1 Sensitivity of the Drift Function 

For forced response problems, the first unsteady quantity that must be examined is 
the rotational velocity, \ R . In Chapter 2 it was shown that the rotational velocity 
is a function of the drift function, A, the stream function, \&, in addition to the 
cascade geometry and prescribed flow conditions. Since the method for computing 
the sensitivity of the cascade geometry and stream function were presented earlier 
in this chapter, only the sensitivity of the drift function needs to be determined to 
calculate the sensitivity of the rotational velocity. 

Since the perturbed grid will still follow the flow streamlines, the sensitivity of the 
drift function may be calculated by linearizing Eq. (3.55), i.e., 


A- • = A- 
*- 


i ,i 


+ 2 


( x i,j ~ a,-w)(g-,j ~ <-u) + (Vi,j ~ Vi-ijWi, j ~ y'i-u) 


$ 






■l,j 


(xjj X;_ij) -)- {yjj yi-ij) 

(** - Si-w) 2 




(4.33) 


Note that the drift function at the inlet boundary remains prescribed to be zero, 
as shown in Eq. (3.55). Earlier in this chapter, the steady flow sensitivity analysis 
procedure was described, which resulted in expressions for the grid perturbations x' 
and y 1 and the steady potential perturbation Hence, the summation in Eq. (4.33) 
may be computed in a straightforward manner. 


4.3.2 Sensitivity of the Rotational Velocity 

The sensitivity of the rotational velocity, v H , may be computed through careful ap- 
plication of the chain rule to the definition of the rotational velocity, Eq. (2.37). 
Considering that the coefficients Ci and c 2 as well as the wave numbers K\ and I< 2 
may have perturbations associated with them in addition to the perturbations al- 
ready discussed in this chapter, the sensitivity of the rotational velocity, v fl/ , may be 
written as 

v R/ = {[(c) + jK[$ + jK \ $') VA + ( Cl +jK 1 $) (V'A + VA') 
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+ 


(^+jK&+jK 2 $' + 


dll' + dH ) 


vt 


4- (c 2 + ]K 2 $ + ||) (V # * + V#') 


x 


+ j {K[ A + K\ A' + K' 2 T + K 2 'f , ) 

(c\ + jK i$) VA + ^c 2 + jK 2 $ + ||) VT | exp[j(I< 1 /L + K 2 1)} (4.34) 


where V 7 refers to the perturbation of the finite difference representation of the gra- 
dient operator due to the grid perturbation. Recall from Chapter 3 that the gradient 
operator used in this report is evaluated using a generalized coordinate transformation 
to the (E, H) grid generation coordinate system [see Eqs. (3.67) and (3.68)]. Hence, 
perturbations in the grid locations result in a perturbation of the discretized gradient 
operator. 


4.3.3 Sensitivity of the Blade Structural Dynamics 

For flutter problems, the first task in the unsteady sensitivity analysis is to determine 
the perturbation of the blade motion. Once the perturbation in the vibratory mode 
shape of the blade is known, the perturbation of the grid motion may be calculated, 
so the complete unsteady sensitivity analysis may be computed. 

Thus far, very little has been said about the structure of the blade. In Chapters 
2 and 3, it was assumed that the vibratory frequency and mode shape are specified 
as part of the flutter problem description. The question then becomes, how do small 
changes in the blade shape affect the frequency and mode shape of the vibration? 

It is beyond the scope of the present analysis to discuss a complete structural 
model of a turbomachine blade. In this analysis it is assumed that there is some 
independent structural model of the blades (e.g., a finite element model) that is 
used to calculate the frequency and mode shape of vibration. For turbomachinery 
blades, unlike aircraft wings, it is reasonable to assume that the ratio of the mass 
of the structure to the mass of the fluid is sufficiently high that the effect of steady 
aerodynamic loading on the calculation of the structural modes may be neglected. 
Hence, the sensitivity of the blade vibration may be calculated independently from 
the present analysis. 

The sensitivity of the equations describing the blade structure may be calculated 
using some of the techniques described earlier or by using similar methods developed 
for analysis of structures [41, 42]. In any event, it is clear that the frequency and 
mode shape of vibration are functions of the design variables. For example, if the 
airfoil is described using the modified NACA four digit definition (see Appendix B), 
the sensitivity of the frequency can be written as 


, _ du du du du du , 

“ “ dm,™ 1 dm, ‘ + dll' + dV‘ + drj 1 


(4.35) 


where m t , rn c , £ t , £ c , and r t are the shape definition variables of the NACA airfoil 
description, as discussed in Appendix B. Similarly, the sensitivity of the mode shape 
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dt ,,df , df X dt , 

d ^, m ‘ + aN' ‘ at/ at /’ + S7* 


To calculate the derivatives in Eqs. (4.35)-(4.36), a sensitivity analysis must be 
performed on the structural model used to compute the nominal natural frequency 
and mode shape of vibration. Once the derivatives of the frequency and mode shape 
have been computed, Eqs. (4.35)-(4.36) may then be substituted into the overall 
unsteady sensitivity analysis expression, Eq. (4.8). 


4.3.4 Sensitivity of the Grid Motion 

Now that the sensitivity of the blade motion has been determined, the next task is 
to compute the perturbation of the grid motion in the computational domain. In 
Chapter 3, the field grid motion was shown to be a linear interpolation of the motion 
on the airfoil surface [Eq. (3.69)]. Hence, the sensitivity analysis is quite simple, 
since the distribution fractions Li and Lj do not change, (i.e., they are not design 
variables). The perturbation of the grid motion, then, may be written as 

f i,i — Li Lj f^iR (4-37) 

where fj^R is the sensitivity of the blade mode shape calculated earlier. 


4.3.5 Sensitivity of the Finite Element Assembly 

Now that all of the perturbation terms have been computed, the right-hand side of 
Eq. (4.8) may be formed. The matrices on the right-hand side of this equation are 
large sparse matrices that may be obtained by linearizing the computer code which 
assembles the global matrix [W] and vector e. The elemental stiffness matrix and 
force vector, along with all of the unsteady boundary conditions, all have a first-order 
sensitivity that may be computed using the methods shown in this chapter. 

For example, the sensitivity of the elemental stiffness matrix, [k] [defined in 
Eq. (3.74)] may be computed by linearizing the matrix with respect to all of the 
variables it is a function of, and multiplying by the actual perturbation, i.e., 




(4.38) 


Since 3?' and x' may be determined from the steady flow sensitivity analysis, and 
u}' may be computed from the structural sensitivity analysis, the sensitivity of the 
stiffness matrix may be computed. To increase the efficiency of the right-hand side 
calculation, the perturbation of the stiffness matrix is multiplied by the appropriate 
entries in the (f) vector at the element level. The details of this procedure do not 
substantially increase the understanding of the sensitivity analysis procedure, and 
are best understood by examining the actual computer code. 

The sensitivity of the force vector, e, [see Eq. (3.75)] may be computed in a similar 



(4.39) 
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Like the sensitivity of the stiffness matrix, the sensitivity of the force vector [i.e., the 
right-hand side of Eq. (4.39)] is best computed at the element level. 


4.3.6 Near-Field Boundary Conditions 

Periodic Boundary Condition 

The periodic boundary condition may contribute inhomogeneous terms to the un- 
steady sensitivity equation because the transformation matrix [I] [see Eq. (3.77)] is 
a function of the interblade phase angle, a. As a result, for example, the linearized 
form of the block equation to be solved at z-stations from i = 2 to i = ILE — 1 [see 
Eq. (3.81)] may be written as 

[ifMPP'i-, + ffl T [b,][i]?. + = (if'e. + [i] T e; 

- ([Ina,][i])X, - (ffl T Ntq)'^. - ([i] T M[i])'5 j+I (4.«) 

Although this appears to be a formidable equation to solve, the sparseness of the [I] 
matrix requires very few entries to be perturbed. Hence, the nominal [a,-], [b t ], and 
[cj] matrices need not be stored in their entirety to perform the sensitivity analysis. 

Airfoil Surface Boundary Condition 

The sensitivity of the unsteady airfoil surface boundary condition need only be con- 
sidered when v R is not zero on the surface. As described in Chapter 3, v R is zero 
on the airfoil surface and wake for airfoils using the Atassi modification of the rota- 
tional velocity. For geometries where this is not the case (e.g., flat plate airfoils using 
the Goldstein formulation), the equation to be solved on the airfoil surfaces may be 
derived by linearizing Eqs. (3.83) and (3.84) about the nominal solution. 


Downstream Wake Boundary Conditions 


Now consider the perturbation of the downstream wake boundary conditions, given by 
Eqs. (3.85)— (3.89) . The wake displacement has a perturbation that must be computed 
along with the perturbation to the unsteady potential. The first boundary condition, 
Eq. (3.85) enforces the pressure continuity across the wake. The appropriate form of 
the perturbation of the pressure continuity equation may be obtained by linearizing 
Eq. (3.85) about the nominal solution. 

The other condition is the conservation of mass on the wake surface, given by 
Eqs. (3.87)-(3.89). In matrix form, the linearized form of the block equation to be 
solved at z-stations from i = ITE + 1 to i = IMAX — 1 [see Eq. (3.91)] may be 
written as 


[a;]0 + [b,]0 . + [cj]0 , +1 = e' - la.]'0,-i - [bi]'0j - [c,]'0, +I 


(4.41) 


Finally, the perturbation of the rotational velocity may be computed in the same 
way as on the airfoil surface, if it is required. 
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4.3.7 Far-Field Boundary Conditions 
Upstream Far-Field Boundary Conditions 

There are new inhomogeneous terms at the far-held boundaries as well. At the up- 


stream far held, the equation to be solved is 

[b\)4>\ + [c,]j/ 2 - ef - [bj]'^ - [c , )'fa (4.42) 

where the inhomogeneous term ef is 

- ([^*] / + z '[^} + 2 [ci]') 4>\ + ([b[] + .z[ci]) 4>\ (4.43) 

and 

[bir^Ibir + Mtl + Ia,]^ (4.44) 

The perturbation of the transition matrix, [T], using the notation of Chapter 3, is 

[T]' = [F 3 ]'[F 2 ][Fi] + [F 3 ][F 2 ]'[Fx] + [F 3 ][F 2 ][F a ]' (4.45) 

~p 

and the perturbation of the particular part of the potential, 4> , is 

(-“[hi] + [bi] + dci]) 4 >\ = 

— (y~ ^[hi] + -[hi]' + [bj]' + z'[c a ] + z[ci]^j (4.46) 


Hence, it is useful to retain the nominal [ax], [Si], and [£x] matrices for the sensitivity 
analysis. 

Downstream Far-Field Boundary Conditions 

Downstream, the equations are somewhat more complicated, since the wake has a 
perturbation associated with it. The equation to be solved at the downstream far- 
held boundary, is 

[a i]4>i-\ + [b /]<£/ = ej' - [a/]V/_x - [b/]'^ (4.47) 

where 

+ Iw + [b;i'j L + (;N + [6;i) (4.48) 

The perturbation of the downstream transition matrix, [T*]', is 

[T*]' = [T]'[I - D] + [T] [I - D]' + z'[D\ + z[ D]' (4.49) 

and the perturbation of the particular part of the potential is 

(—“[a/] + [b /] + z[ c/]) 4> ! = e'j 

- ^[a/] + “[a/] 7 + [b /]' + z'[ c/] + z[ c/]'j fa (4.50) 

The nominal [a./], [S /], and [c/] matrices are also stored for use in the sensitivity 
analysis. 


Chapter 5 
Results 


This chapter will illustrate how the present method may be used to analyze typical 
aeroelastic and aeroacoustic problems encountered in turbomachines. The results 
of the sensitivity analysis will be used to suggest design changes to improve the 
unsteady aerodynamic behavior. The newly redesigned airfoils will then be analyzed 
and compared to the initial design to assess the effectiveness of the procedure. The 
computational cost of the procedure will also be examined. In Section 5.1, the flutter 
stability of a compressor rotor will be analyzed as an example of a typical aeroelastic 
problem. The aeroacoustic capabilities of the present method will be demonstrated in 
Sections 5.2 and 5.3, where the acoustic response of two fan exit guide vanes (EGV) 
due to an upstream vortical rotor wake will be examined. Finally, Section 5.4 will 
summarize the conclusions from the results of the analysis. 

In the computational results presented here, lengths are nondimensionalized by the 
airfoil chord, c, steady and unsteady velocities by the upstream steady velocity |V_oo|, 
frequencies by IV-ool/c, pressures by .ft-oolVioJ, lifts per unit span by R-ooIVloJc, 
and moments per unit span by ALoojVioJc 2 . 

5.1 Aeroelastic Analysis and Design of a Com- 
pressor 

5.1.1 Steady Flow Through a Compressor 

To demonstrate the aeroelastic capabilities of the present sensitivity analysis, a linear 
cascade of NACA four digit series airfoils will be analyzed. The nominal cascade 
considered here is similar to modern compressor cascades, and is composed of NACA 
5506 airfoils. The inflow Mach number M_ <*, is 0.5, the inflow angle Cl-oo (measured 
from the axial direction) is 55°, the stagger angle 0 is 45°, and the blade-to-blade 
gap G is 0.9. 

The steady flow through the cascade was computed using two different computa- 
tional grids — a 65 x 17 node H - grid and a 129 x 33 node grid. Figure 5.1 shows the 
nominal steady surface pressure, P, for these two different grid resolutions. Note the 
good agreement between the coarse grid and fine grid solutions, indicating that the 
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Figure 5.1: Steady surface pressure of cascade of NACA 5506 airfoils. M _ TO = 0.5, 
fi.oo = 55°, 0 = 45°. 

solution is “grid converged.” The steady surface pressure may be integrated along the 
chord to determine the steady forces acting on the airfoil. In this case, the steady lift, 
L (measured normal to the airfoil chord), is 0.2907, and the steady drag, D (measured 
tangent to the airfoil chord), is —0.0177. The steady moment may be calculated by 
integrating the product of the pressure and the distance to the point about which the 
moment is being taken. For this case, the moment measured about the leading edge, 
Mle, is —0.1215. It should also be noted that the flow is entirely subsonic with a 
maximum Mach number on the suction surface of about 0.61. 

5.1.2 Unsteady Flow Through a Compressor 

Now that the steady flow through the cascade has been computed, we consider the 
unsteady flow due to plunging and torsional vibration of the airfoils. On the surface 
of the reference (i.e., unrotated) airfoil, the grid motion vector f for this rigid body 
motion may be written as 

f (x,y) = [-(y - y„)ar]i + [(x - x p )a + h] j (5.1) 

where h is the amplitude of the bending vibration normal to the blade chord, a is 
the amplitude of the pitching motion (positive nose down), and the airfoil is pitching 
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about the point ( x p , y p ). The vectors i and j are unit vectors in the x- and y-directions, 
respectively. 

As described in Chapter 1, the unsteady aerodynamic portion of a flutter analysis 
is to determine the unsteady work done on the blade due to self-induced oscillations. 
The blade work is the product of the aerodynamic force and the unsteady blade 
motion. A particularly useful quantity is the unsteady work done on the blade over 
one period of the oscillation, denoted Wcycle- The complete work per cycle is defined 
as [55] 

/■ A r - di A 

Wcycle = ~ j q f r ' n ds dt (5.2) 

where A is the period of the oscillation, P is the complete time- varying pressure, and 
n is the unsteady unit normal to the blade surface. A positive value of the work 
per cycle indicates that the fluid is adding work to the motion, which continues to 
increase the amplitude of the motion, and is therefore considered to be destabilizing. 
After some algebra [56], it may be shown that 

wcycle = * f [(pK f / ~ P/ f «) • n - P (f R ■ n'/ - fj • n'*)] ds (5.3) 

where the subscripts R and I refer to the real and imaginary parts of the quantity, 
respectively. The vector n' arises from the tilting of the unit normal due to the 
pitching of the blade. On the reference airfoil, note that the grid motion is defined 
to be purely real. We will define the unsteady lift, I, to be defined as 

i — —7 xh <J> p dx — a j) P dx (5-4) 

In a similar fashion, the unsteady moment, mx, will be defined as 

m T = 7r a jf p [(x - x v )dy + (y - y p )dx } (5.5) 

We wish to express the work per cycle in terms of the unsteady lift and moment. 
To accomplish this, consider a single-degree-of-freedom bending or torsional vibration 
(i.e. , the blade is either pitching or plunging, but not both). In this case, from 
Eq. (5.3) it may be shown that for a purely bending vibration, 

wcycle = tt/i Im[£] (5-6) 

Similarly, for a purely torsional vibration, 

wcycle = *oi I m[m T ] (5.7) 

So for a single degree of freedom motion, it is the out-of-phase (with blade displace- 
ment) component of the lift for plunging motion, or moment for a pitching motion 
that determines whether the motion is stable or unstable. In structural dynamic 
terms, Carta [57] has shown that for single-degree-of-freedom motions the work per 
cycle may be expressed as the aerodynamic damping of the system. Since the struc- 
tural damping of the aeroelastic system is usually small, the aerodynamic damping 
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Figure 5.2: Aerodynamic damping of cascade of NACA 5506 airfoils vibrating in 
plunge at frequencies of 0.4, 0.8, and 1.6 for a range of interblade phase angles. 

essentially determines the aeroelastic stability of the system. For a purely bending 
vibration, the aerodynamic damping H b may be written as 

Hb = t^cycle = — Im[£] 

7 r/i 

and for a purely torsional vibration, Et is 

= w cycle = ~Im[mr] 

7T a 

These expressions clearly show that positive aerodynamic damping indicates that the 
motion is stable. 

To evaluate the stability of the cascade under consideration, we will calculate the 
aerodynamic damping for both plunging and pitching for three different “reduced” 
frequencies. The term “reduced frequency” indicates that the frequency of vibration 
has been nondimensionalized by the airfoil chord, c, and the upstream freestream 
velocity, |V-oo|- Figure 5.2 shows the aerodynamic damping Eg of the cascade vi- 
brating in plunge at three reduced frequencies and for a range of interblade phase 
angles. Note that for plunging motion, the system is stable, that is, the aerody- 
namic damping is positive for all interblade phase angles. However, the aerodynamic 
damping is generally less for low reduced frequencies. The pronounced peaks in the 
damping curves correspond to acoustic resonances of the duct containing the blade 
row. At interblade phase angles between these acoustic resonance points, at least 


(5.9) 


(5.8) 
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Interblade Phase Angle (deg) 

Figure 5.3: Aerodynamic damping of cascade of NACA 5506 airfoils pitching about 
their midchords at frequencies of 0.4, 0.8, and 1.6 for a range of interblade phase 
angles. 

one of the pressure modes propagates unattenuated, i.e., it is “cut-on.” At interblade 
phase angles outside of this region, all the pressure modes decay as the propagate, or 
are “cut-off.” 

Figure 5.3 shows the aerodynamic damping Et for the case where the airfoils 
vibrate in pitch about their midchords. Again, the cascade is least stable at the low 
reduced frequencies. In particular, note that the system is unstable (E t < 0) for 
several interblade phase angles at the lowest reduced frequency u> of 0.4. In other 
words, we have discovered an instability in the aeroelastic system. 

We wish to use the sensitivity analysis developed in this report to suggest changes 
in the shape of the airfoil so that the instability may be eliminated. Hence, it is 
useful to examine in detail the case where the motion is least stable. As shown in 
Fig. 5.3, the least stable case occurs when the airfoils pitch about their midchords 
with a reduced frequency u of 0.4 and an interblade phase angle a of 60° (this is the 
least stable interblade phase angle for the reduced frequency u of 0.4). 

Figure 5.4 shows the real and imaginary parts of the complex amplitude of the 
nominal unsteady pressure, p, on the surface of the reference airfoil computed using 
two different grid resolutions, a 65 x 17 node grid and a 129 x 33 node grid. Note 
that the imaginary part of the pressure distribution is the part that does work on the 
vibrating airfoil. For this case, it is apparent that the imaginary part of the pressure 
difference across the airfoil is generally negative over the front half of the airfoil and 
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Distance Along Chord, ^/c 


Figure 5.4: Imaginary part of unsteady surface pressure of NACA 5506 airfoils pitch- 
ing about their midchords. u> = 0.4, a = 60°. 
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positive over the aft half. Thus, since the airfoil pitches about its midchord (positive 
nose up), the unsteady pressure does positive aerodynamic work (corresponding to 
negative aerodynamic damping) on the airfoil over most of the airfoil, making the 
cascade unstable in pitch for these conditions. 

5.1.3 Sensitivity Analysis 

Having computed the nominal flow through the cascade, the effect of five different 
design parameters on the steady flowfield will now be studied. Two of these pa- 
rameters are from the NACA four digit airfoil definition, i.e., the magnitudes of the 
airfoil thickness and camber. (See Appendix B for a discussion of the NACA airfoil 
definition.) Each of these quantities are measured in fractions of the airfoil chord c. 
Also considered are the effect of changes in the cascade stagger angle 0, and blade- 
to-blade gap G. Finally, an additional design variable, the reflex, is introduced here. 
The reflex parameter modifies the height of the mean line 

h m = h c + rcsin(27r£/c) (5.10) 

where h m is the height of the mean line, h c is the height of the mean line due to 
camber, r is the magnitude of the reflex in fractions of chord, c is the chord, and is 
the distance along the airfoil chord. Reflex has been added to provide an “S shape” 
parameter to the design, to allow more flexibility in the airfoil shape than the strict 
NACA definition. 

Figure 5.5 shows the sensitivity of the steady surface pressure to changes in these 
five geometry variables. The sensitivities are computed using the present sensitivity 
analysis; all results were computed on a 65 x 17 node grid. To check these results, the 
sensitivities using a finite difference approach are also computed. The finite- difference 
result is computed by solving for the steady flow about two slightly different airfoils, 
differencing the two solutions, and dividing the result by the difference in the airfoil 
parameter. Note the excellent agreement between the two solutions indicating that 
the effect of small changes in the design variables is linear, and that the present 
sensitivity analysis correctly predicts the sensitivities. Also, not surprisingly, the 
largest sensitivity in pressure occurs near the leading edge of the airfoil. 

Next, the surface pressure sensitivities were integrated to obtain the sensitivity of 
the steady lift and drag (measured normal to and along the chord) and the moment 
about the leading edge. These results are given in Table 5.1. The sensitivity to 
changes in maximum camber location is also shown. In addition, the sensitivity of 
the lift in the y-direction (the cascade direction) is tabulated. The steady lift in the 
y-direction is a measure of the turning done by the cascade and hence is related to the 
steady work done by the cascade. Table 5.1 shows that the lift in the y-direction is 
most sensitive to changes in camber, stagger angle, and reflex. Since these parameters 
control the metal angle of the trailing edge, and the deviation between the exit flow 
angle and the metal angle is small for cascades, one would expect them to have a 
strong influence on the steady lift. 

Next, the sensitivities of the unsteady surface pressure to changes in geometry are 
computed. Figures 5.6 and 5.7 show the real and imaginary parts, respectively, of 


C-Z, 
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Figure 5.5: Sensitivity of steady surface pressure of cascade of NACA 5506 airfoils to 
perturbations in thickness, camber, stagger, gap, and reflex. M_ = 0.5, O-oo = 55°, 
0 = 45°. 
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Table 5.1: Sensitivity of steady forces and moment. The nominal steady lift, L, is 
0.2907, the nominal drag, D, is —0.0177, the nominal moment about the leading edge, 
Mle , is —0.1215, and the nominal lift in the j/-direction, Ly, is 0.1931. 


Design Variable 

L' 

D' 

M' le 

L'y 

Thickness 

-0.1935 

-0.0135 

-0.1234 

-0.1464 

Camber 

1.5637 

0.1446 

-1.4003 

1.2080 

Stagger 

-0.6632 

-0.0642 

-0.0548 

-0.5144 

Gap 

0.2743 

-0.0244 

-0.0764 

0.1767 

Max Camber Location 

0.0873 

0.0083 

-0.1060 

0.0676 

Reflex 

-1.5506 

-0.1476 

1.9235 

-1.2008 


the sensitivity of the unsteady pressure to small changes in six design variables (the 
reduced frequency u is included as a design variable for unsteady flow calculations). 
All results were computed on a 65 x 17 node grid. The sensitivities are also compared 
to a finite difference calculation. Note the excellent agreement between the two so- 
lutions indicating that the present method correctly predicts the sensitivities. Also, 
the imaginary parts of the sensitivities to changes in stagger and reflex have pressure 
distributions that are fairly large in magnitude and have shapes that would tend to 
do work on pitching airfoils. That is, the sign of the pressure difference across the 
airfoil changes at roughly the midchord of the airfoil. 

The sensitivities of the surface pressure to design variables may now be integrated 
to obtain the sensitivities of the aerodynamic damping. Table 5.2 shows the sensitiv- 
ity of the aerodynamic damping to small changes in the design variables, including 
maximum camber location. The column labeled “Unconstrained” gives the sensitivity 
of the aerodynamic damping to changes in a single parameter. Here 5 'y is the sensi- 
tivity of the aerodynamic damping due to pitching motions, and E ' b is the sensitivity 
of the aerodynamic damping due to plunging motion. In both cases, the nominal 
reduced frequency u is 0.4. Note that as expected, stagger and reflex have a strong 
influence on the aerodynamic damping in pitch. Also note that for both pitching and 
plunging, the sensitivity of the damping to changes in frequency is positive. This is 
consistent with the results shown in Figures 5.2 and 5.3. 

The results in the “Unconstrained” column of Table 5.2, however, can be somewhat 
misleading since changing each design variable independently also changes the steady 
work done by the blade row and changes the steady incidence at the leading edge 
of the airfoils. Generally, one would want to leave these quantities unchanged. To 
avoid this difficulty, it is useful to let two of the design variables “float” so that the 
steady turning done by the blade row, L y /G, and the leading edge incidence angle, 
a, remain constant. In this study, the stagger angle 0 and the reflex r are allowed to 
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Figure 5.6: Real part of sensitivity of unsteady surface pressure of NACA 5506 airfoils 
pitching about their midchords due to perturbations in thickness, camber, stagger, 
gap, reflex, and frequency, u = 0.4, a = 60°. 
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Distance Along Chord.^/c 


Figure 5.7: Imaginary part of sensitivity of unsteady surface pressure of NACA 5506 
airfoils pitching about their midchords due to perturbations in thickness, camber, 
stagger, gap, reflex, and frequency, ui = 0.4, a = 60°. 
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Table 5.2: Sensitivity of aerodynamic damping. The nominal aerodynamic damping 
in torsion, Zt, is —0.0214, and the damping in plunging, Zb, is 0.8882. 



Unconstrained 

Constrained 

Design Variable 


— B 

“T 

'-‘B 

Thickness 

-0.2208 

0.3046 

0.1344 

-0.1879 

Camber 

-0.0018 

-1.9868 

-5.2561 

3.3668 

Stagger 

-0.6672 

1.4079 

— 

— 

Gap 

0.1723 

-0.0237 

0.2643 

-0.1511 

Max Camber Location 

-0.0383 

-0.1143 

0.0300 

-0.0159 

Reflex 

0.7362 

2.0131 

— 

— 

Frequency 

0.4030 

1.3337 

0.4030 

1.3337 


float. For example, then, if the gap G varies, the stagger angle and reflex must vary 


such that 

O' + r' + K? G‘ L *ff = 0 

30 dr oG G 

(5,11) 

and 

da a> GL da P± da r’ -n 

se @ + a7 r + sg g "° 

(5,12) 


Equations (5.11) and (5.12) give two equations for the two unknowns 0' and r' in 
terms of the gap perturbation G' and the sensitivities. The sensitivity of the incidence 
angle, a , and the resulting perturbations in 0 and r for each design variable are shown 
in Table 5.3. In Table 5.2, the column labeled “Constrained” refers to the sensitivities 
to each variable using this procedure. For both the pitching and plunging cases, it is 
clear that changing the camber has a very strong effect on the aerodynamic damping. 
In the pitching case, an increase in camber is destabilizing; in the plunging case, an 
increase in camber is stabilizing. 

5.1.4 Redesign of a Compressor for Aeroelastic Stability 

Next, the constrained sensitivity analysis is used to redesign an unstable cascade 
to make it stable. The nominal cascade has a reduced frequency u of 0.4 and an 
interblade phase angle a of 60°. Note from Table 5.2 that decreasing the camber has 
a stabilizing influence on torsional flutter. Thus, for the first redesign (Redesign A), 
the camber is reduced by 0.004 units. Using the constraint relations shown earlier, 
this requires that the stagger angle must be reduced by approximately 1.4° and the 
reflex must be increased by 0.0064 units (see Table 5.3). Although the sensitivity 
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Table 5.3: Sensitivity of incidence angle to design variables, and resulting perturba- 
tions in stagger angle and reflex for a unit change in design variables. 


Design Variable 

da/d(v ar) 

0' 

r' 

Thickness 

0 

-0.4533 

0.0721 

Camber 

2/4 

6.1033 

-1.6080 

Stagger 

1 



Gap 

0 

0.5465 

-0.0870 

Max Camber Location 

— 2m c /(4) 2 

-0.0274 

0.0680 

Reflex 

2t r 



Frequency 

0 

0.0000 

0.0000 


analysis predicts that these changes alone will make the airfoil stable, the sensitivity 
analysis also predicts a large steady pressure gradient on the suction surface near 
the leading edge which would very likely cause the flow to separate, an undesirable 
result. To reduce the adverse pressure gradient in the pressure distribution, the 
thickness is increased by 0.02 units, which in turn requires the stagger be reduced by 
approximately 0.52° and 0.0014 units of reflex be added. 

For the second redesign (Redesign B), the gap G is increased by 0.1. Again 
Table 5.2 predicts that this change will make the cascade stable, and requires that 
the stagger angle be reduced by approximately 0.61° and 0.0017 units of reflex be 
added. 

Figure 5.8 shows the computed steady surface pressure on the nominal and re- 
designed airfoils. Also shown is the pressure predicted by the linear sensitivity analy- 
sis. The good agreement between the two indicates that steady nonlinear geometrical 
effects are small, at least for the subsonic flow conditions considered here. For Re- 
design A, although the steady lift on the airfoil in the y-direction has only slightly 
changed, the pressure distribution has changed significantly. For Redesign B, the net 
lift in the y-direction has increased, since the increase in gap means that the steady 
work per airfoil must increase. Note that the pressure gradient on the suction surface 
is larger for both redesigned airfoils. Both redesigns are therefore likely to increase 
somewhat the aerodynamic losses of the cascade. 

Figures 5.9 and 5.10 show the real and imaginary parts of the unsteady pressure 
on the surface of the redesigned airfoils for Redesign A and Redesign B, respectively. 
Although in both cases the real part shows very little change, the imaginary part 
shows significant changes, particularly on the suction surface. Although there is a 
larger difference between the sensitivity analysis prediction and the actual pressure 
distribution than in the steady case, the sensitivity analysis still provides an excellent 



102 


CHAPTERS. RESULTS 



CM 



Figure 5.8: Steady surface pressure of cascade of redesigned airfoils. M- <*, = 0.5, 
fl_oo = 55°. 

qualitative estimate of the actual unsteady flow behavior. The actual computed 
damping of the Redesign A cascade is 0.0086, indicating that the new cascade is 
stable. The damping of the Redesign B cascade is 0.0015, so this cascade is also 
stable. 

At this point, it is useful to consider the accuracy of the sensitivity analysis pre- 
diction for relatively large changes in the airfoil shape. If the change in the design 
variables in Redesign A is considered to be one unit, Fig. 5.11 shows the change in the 
aerodynamic damping predicted by the sensitivity analysis and the actual change for 
various magnitudes of the design change. Figure 5.11 shows that the behavior of the 
aerodynamic damping is remarkably linear, i.e., the airfoils become physically unre- 
alistic before substantial differences occur between the sensitivity analysis prediction 
and the actual airfoil behavior. 

Figures 5.12 and 5.13 show the aerodynamic damping of the redesigned airfoils 
for a reduced frequency u of 0.4 for a range of interblade phase angles a . In each 
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Redesign A 
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Figure 5.9: Real and imaginary parts of unsteady surface pressure of redesigned 
airfoils (Redesign A) pitching about their midchords. u> = 0.4, a = 60°. 
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Figure 5.10: Real and imaginary parts of unsteady surface pressure of redesigned 
airfoils (Redesign B) pitching about their midchords. u> = 0.4, a = 60°. 
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Figure 5.11: Accuracy of sensitivity analysis for various perturbation amplitudes. 


figure, the original nominal damping, the damping of the redesigned airfoils predicted 
by the sensitivity analysis, and the actual damping of the redesigned airfoils are 
presented. Note that both redesigned airfoils are stable for all interblade phase angles. 
In addition, the sensitivity analysis prediction gives excellent estimates of the actual 
damping of the redesigned airfoils. 


5.1.5 Computational Efficiency 

Finally, a note about computational times. Table 5.4 shows the CPU time required to 
perform various calculations using the present method on a Silicon Graphics Indigo 
R4400 workstation. All time calculations were performed using a 129 x 49 node 
computational grid. The steady sensitivity analysis requires only a fraction of the 
CPU time necessary to perform a single nominal steady calculation. For the six 
design variables considered here, the unsteady sensitivity analysis required about 
three times the CPU time as a single nominal unsteady calculation, but only about 
one-fifth of what was required for a finite difference sensitivity analysis. Furthermore, 
the present sensitivity analysis, unlike the finite difference analysis, is not susceptible 
to truncation and round-off errors. 

It should be noted, however, that for this aeroelastic example, the computational 
grid was computed independently from the steady potential, since a streamline grid is 
not necessary for this case. Decoupling the grid and steady flow equations significantly 
reduces the computational expense of the present method, as will be shown in the 
next section. 
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Figure 5.12: Aerodynamic damping of cascade of redesigned airfoils (Redesign A) 
pitching about their midchords at a frequency of 0.4 for a range of interblade phase 
angles. 


Table 5.4: Computational times for present method using 129 x 49 node grid. 


Procedure 

CPU Time (sec) 

Nominal Steady 

85.3 

Nominal Unsteady 

12.1 

Steady Sensitivity Analysis (5 var) 

10.2 

Unsteady Sensitivity Analysis (6 var) 

30.6 

Finite Difference Steady Sensitivity Analysis (5 var) 

853.0 

Finite Difference Unsteady Sensitivity Analysis (6 var) 

145.2 
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Figure 5.13: Aerodynamic damping of cascade of redesigned airfoils (Redesign B) 
pitching about their midchords at a frequency of 0.4 for a range of interblade phase 
angles. 


5.2 Aeroacoustic Analysis and Design of a Fan 
Exit Guide Vane 

To demonstrate the aeroacoustic capability of the present method, in this section a 
cascade of exit guide vanes (EGV) typical of those found in modern high-bypass ratio 
fans is analyzed. The nominal airfoil shape is a NACA 8508-65 profile. The ratio of 
the number of fan rotor blades to EGVs, Nr/Nv , is 0.4. The blade-to-blade gap G 
is 1.0, the inlet Mach number M - «, is 0.5, the inlet flow angle is 30°, and the 
stagger angle 0 is 16°. The wheel speed of the upstream rotor VR otor is 1.5. 

The physical system being modeled here is shown schematically in Fig. 5.14. A 
row of exit guide vanes (EGVs) is subjected to unsteady aerodynamic excitation 
arising from interaction with the viscous wakes from the upstream rotor. In the 
rotor (relative) frame of reference, these wakes are steady, and the steady freestream 
velocity is Vp. e i. However, in the EGV (absolute) frame of reference, the wakes appear 
to be unsteady, and the steady freestream velocity is V^bs- The two reference frames 
are related by the wheel speed, VR 0 ton due to the relative motion of the rotor. 

The unsteadiness of the rotor wakes may be decomposed into harmonics with 
frequencies that are multiples of blade passing frequency (BPF). The BPF is referred 
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Figure 5.14: Schematic showing wake/EGV interaction. Inserts in blade passage show 
contours of “drift” (top) and a streamline computational grid (bottom). 

to as the fundamental (i.e., 0th harmonic) frequency that the rotor wakes impinge 
on the stator. Typically, the number of rotor blades and stator vanes are chosen so 
that no pressure waves propagate (“cut-on”) at this frequency. Unfortunately, it is 
not possible to choose the blade counts so that pressure waves do not propagate at 
harmonics of the BPF. Hence, the goal for this investigation is to use the sensitivity 
analysis to suggest design changes in the airfoil shape to reduce the magnitude of the 
propagating pressure waves at multiples of BPF. 

5.2.1 Steady Flow Through a Fan Exit Guide Vane 

First, consider the steady flow through the EGVs. The steady flow was computed 
using an H - grid containing 129 nodes in the streamwise direction and 49 nodes in 
the normal direction (a 129x49 grid). Figure 5.15 shows the nominal steady surface 


5.2. AEROACOUSTIC ANALYSIS AND DESIGN OF A FAN EGV 


109 



Symbols 

Present Method 

Lines 

Linearized Euler 


CL 



0.0 0.2 0.4 0.6 0.8 1.0 

Distance Along Chord, x/c 


Figure 5.15: Steady surface pressure for cascade of NACA 8508-65 airfoils. 0 = 16°, 
O.oo = 30°, M_oo = 0.5, G = 1.0. 

pressure P for this case. The steady lift, L , is 0.3942, and the steady drag, D, 
is —0.0168. The moment measured about the leading edge, Mle, is —0.1837. For 
comparison, Fig. 5.15 also shows the grid-converged solution computed using an Euler 
code [30]. The excellent agreement between the two solutions, while reassuring, is to 
be expected since the steady flow is subsonic, irrotational, and homentropic. 

5.2.2 Unsteady Flow Through a Fan Exit Guide Vane 

Next, consider the acoustic response due to viscous wakes from the upstream fan im- 
pinging on the EGV. When viewed in the EGV frame of reference, the wake excitation 
has temporal frequencies u> n = 27rnVR.ot.or/GK, where VR otor is the wheel speed of the 
fan, Gr is the blade-to-blade gap of the rotor, and n takes on all integer values. Blade 
passing frequency (BPF) corresponds to n = 1. The corresponding interblade phase 
angles are cr n = — 27T hG/Gr. The acoustic response of this cascade to excitations 
at lx BPF is cut-off, i.e. the unsteady pressure decays exponentially away from the 
EGVs. For excitations at 2xBPF, however, a single acoustic mode with interblade 
phase angle a = 72° is cut-on in the upstream and downstream regions. 

In the following, the acoustic response of the EGV to 1 xBPF and 2xBPF vortical 
gusts with “unit amplitude” are considered. A unit amplitude gust is one in which the 
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magnitude of the perturbation velocity normal to the steady flow direction would be 
unity at the leading edge of the airfoil if the steady flow were uniform and undeflected 
by the EG V. The real and imaginary parts of the unsteady surface pressure computed 
using the present analysis are shown in Figs. 5.16 and 5.17. Also shown for compari- 
son is the pressure distribution computed using the linearized Euler analysis of Hall 
and Clark [30]. The agreement between these two theories is quite good, especially 
considering the high reduced frequency (co = 7.54) of the 2xBPF case. These results 
indicate that the 129x49 grid is sufficiently fine to resolve the acoustic behavior of 
the EGV — at least up to 2xBPF. They also indicate that the rotational velocity 
has been formulated and computed correctly. 

Figures 5.18 and 5.19 show contours of the magnitude of the unsteady pressure 
for the lxBPF and 2xBPF cases, respectively. In both cases, note that the unsteady 
pressure is somewhat larger in the downstream region than in the upstream region. 
The unsteady pressure in the far field of the 2xBPF case was Fourier transformed in 
the circumferential direction to determine the magnitude of the single cut-on pressure 
wave upstream and downstream (cr = 72°). For this case, the magnitude of the cut- 
on unsteady pressure wave (for a unit strength gust) is 0.102 upstream and 0.249 
downstream. 

5.2.3 Sensitivity Analysis 

Next, the change in the steady and unsteady aerodynamic response due to small 
changes in the EGV geometry is considered. The steady sensitivity analysis was 
performed using eight different design variables. Five of these variables correspond 
to the modified NACA definition parameters: magnitude of thickness, maximum 
thickness location, magnitude of camber, maximum camber location, and leading 
edge radius. Two design variables define the placement of the airfoils in the cascade: 
the blade-to-blade gap G and stagger angle 0. The final design variable is the reflex, 
which was defined earlier. 

Using the steady sensitivity analysis outlined in the previous section, the sensitiv- 
ity of the lift, drag, and moment to changes in the design variables were computed. 
These results are given in Table 5.5. Like the previous example, Table 5.5 shows 
that the steady work is very sensitive to changes in camber, reflex, and stagger angle. 
This is not too surprising since these parameters affect the trailing edge metal angle 
of the airfoil and hence the turning (work). On the other hand, the steady work is 
insensitive to changes in the maximum thickness location and the size of the leading 
edge radius. 

Next, the sensitivities of the unsteady pressure to small changes in the design 
variables for the 2xBPF case were computed. In addition to the eight design vari- 
ables described above, the frequency w and interblade phase angle cr of the excitation 
have been added. Figures 5.20 and 5.21 show the real and imaginary parts, respec- 
tively, of the sensitivity of the unsteady surface pressure to changes in five of these 
design variables. Also shown are the sensitivities calculated using a finite-difference 
approach, i.e. the difference of the pressures calculated using the nominal analysis on 
two slightly different airfoils normalized by the difference in geometry. The nearly 
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Figure 5.16: Unsteady surface pressure of cascade of NACA 8508-65 airfoils due to 
incoming vortical gust at lxBPF. u = 3.7687, a = —144°. 
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Figure 5.17: Unsteady surface pressure for cascade of NACA 8508-65 airfoils due to 
incoming vortical gust at 2xBPF. u> = 7.54, cr = —288°. 
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Figure 5.18: Contours of magnitude of unsteady pressure for cascade of NACA 8508- 
65 airfoils due to incoming vortical gust at lxBPF. u> = 3.77, a = —144°. 
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gure 5.19: Contours of magnitude of unsteady pressure for cascade of NACA 8508- 
airfoils due to incoming vortical gust at 2xBPF. u = 7.54, a = —288°. 
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Table 5.5: Change in steady flow quantities due to unit perturbations in ten design 
variables. The nominal steady lift, L , is 0.3942, the nominal drag, D, is —0.0168, the 
nominal moment about the leading edge, Mle , is —0.1837, and the nominal lift in 
the y-direction, Ly , is 0.3743. 


Design Variable 

Ly 

L' 

D' 

M' le 

Thickness 

0.0289 

0.0284 

0.0059 

0.1751 

Max Thickness Location 

0.0084 

0.0082 

0.0017 

-0.0164 

Camber 

1.9254 

1.8941 

0.3799 

-1.6300 

Max Camber Location 

0.1991 

0.1959 

0.0394 

-0.2154 

Leading Edge Radius 

0.0000 

0.0000 

0.0000 

-0.0001 

Stagger 

-0.7558 

-0.7266 

0.2448 

0.0392 

Gap 

0.3008 

0.3220 

-0.0313 

-0.1123 

Reflex 

-2.2190 

-2.1827 

-0.4387 

2.4383 


exact agreement between the two methods demonstrates that the present sensitivity 
analysis has been formulated correctly. 

Of particular interest in aeroacoustic applications is the influence of the blade 
shape on the sound radiated upstream and downstream of the EGV. The sensitivity 
of the magnitude of the cut-on propagating pressure waves to changes in the design 
variables is given in Table 5.6. The label “Unconstrained” indicates that each of the 
design variables is perturbed independently. Like the previous example, to compute 
the “Constrained” sensitivity, each of the design variables is perturbed as before, 
but the stagger angle and reflex are allowed to float to satisfy the constraints that 
the steady work and incidence angles remain fixed. Examining the “Constrained” 
columns in Table 5.6, it is clear for example that moving the position of the maximum 
camber location aft now increases the magnitude of the upstream pressure wave and 
decreases the magnitude of the downstream pressure wave. Note further that both 
the gap G and camber have a strong influence on the strength of the acoustic waves 
(although large changes in the gap are more realistic than large changes in camber). 

5.2.4 Redesign of an EGV for Reduced Acoustic Response 

Next, the constrained sensitivities were used to guide the redesign of the EGV to 
reduce the sound pressure levels in the downstream region. As previously noted, the 
camber and gap strongly influence the outgoing pressure waves. Therefore, to reduce 
the downstream pressure wave, the blade-to-blade gap is increased by 0.1 and the 
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Figure 5.20: Sensitivity of unsteady surface pressure on NACA 8508-65 airfoils due 
to perturbations in thickness, camber, stagger, gap, and reflex, u = 7.54, a = —288°. 
□ , suction surface; A, pressure surface. 
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Figure 5.21: Sensitivity of unsteady surface pressure on NACA 8508-65 airfoils due 
to perturbations in thickness, camber, stagger, gap, and reflex, u = 7.54, a = —288°. 
□ , suction surface; A, pressure surface. 
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Table 5.6: Change in unsteady flow quantities due to unit perturbations in ten design 
variables. The nominal magnitude of the upstream pressure wave, |p up |, is 0.102, and 
the magnitude of the downstream pressure wave, |pa own |, is 0.249. 



Unconstrained 

Constrained 

Design Variable 

|Pup|' 

|Pdown | 

bu P r 

Ipdown | 

Thickness 

-0.8985 

0.2494 

-0.9447 

0.3099 

Max Thickness Location 

0.0084 

0.0082 

-0.0810 

0.1161 

Camber 

0.4441 

-1.4103 

-6.0899 

6.5867 

Max Camber Location 

-0.1312 

0.1224 

0.1023 

-0.0953 

Leading Edge Radius 

0.0023 

0.0023 

0.0023 

0.0023 

Stagger 

-0.3471 

0.5910 

— 

— 

Gap 

1.3244 

-1.6624 

1.4421 

-1.8163 

Reflex 

1.8692 

-1.5834 

— 

— 

Frequency 

0.1696 

-0.2830 

0.1696 

-0.2830 

Interblade Phase Angle 

-0.0822 

0.0222 

-0.0822 

0.0222 


camber is decreased by 0.0025. To satisfy the constraints, the stagger angle, 0, must 
be reduced by 2.23°, and 0.0078 units of reflex must be added. Also note that the 
interblade phase angle, cr, is decreased by 0.5027 radians since the interblade phase 
angle is proportional to the gap. For the 2xBPF case, this means that the interblade 
phase angle is decreased from —288° to —317°. 

Figure 5.22 shows the nominal and redesigned airfoils and their respective steady 
surface pressure distributions. The redesigned pressure distribution was computed in 
two ways: first by linear extrapolation using the sensitivities, and second using the 
nonlinear steady flow analysis with the actual redesigned airfoil geometry. The good 
agreement between the two techniques shows that the present sensitivity analysis is 
valid for moderate changes in geometry. Also, one can see in Fig. 5.22 that the net 
lift has increased to account for the additional steady work done by each blade due 
to the increase in blade-to-blade gap. 

Figure 5.23 shows the real and imaginary parts of the unsteady surface pressure 
for the nominal and redesigned cascade. Note that although the sensitivity analysis 
predicts the trends of the actual new unsteady pressure distribution, there is a signif- 
icant difference in magnitude between the two curves. The reason for this difference 
can be explained by examining the size of the design change. Figure 5.24 shows the 
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Figure 5.22: Steady surface pressure for cascade of redesigned airfoils. M_oo = 0.5, 
ft-oo = 30°. 

sensitivity analysis prediction of the change in the upstream and downstream pres- 
sure waves versus the actual change for several different design change amplitudes. 
The design change amplitudes are normalized by the amplitude described above, i.e., 
the above design change is considered to have an amplitude of unity. Figure 5.24 
shows that the differences between the predicted pressure distribution and the actual 
distribution are due to nonlinear effects that are not predicted using the present anal- 
ysis. If a smaller change in the geometry had been made, the prediction would have 
been more accurate. Even for this large a design change, however, the sensitivity 
analysis still gives an excellent qualitative prediction of the change in the magnitude 
of propagating pressure waves. 

Figures 5.25 and 5.26 show the unsteady pressure contours for the redesigned 
cascade for both cases. Comparing with Figs. 5.18 and 5.19, note that the downstream 
pressure levels have been substantially reduced, although the upstream pressure levels 
have been increased. For the 2xBPF case, the magnitude of the downstream pressure 
wave was reduced about 60%, from 0.249 to 0.099. In acoustic terms, the sound 
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Figure 5.23: Real and imaginary parts of unsteady surface pressure for redesigned 
airfoils due to incident vortical gust at 2xBPF. u = 7.54, a = —317°. 





5.2. AEROACOUSTIC ANALYSIS AND DESIGN OF A FAN EGV 


121 



Figure 5.24: 


Accuracy of sensitivity analysis for various perturbation amplitudes. 


pressure level (SPL) is reduced by 8.0 dB, since the change in SPL may be defined as 


Change in SPL = 201og 10 


| Predesign 1 
I Pnomin.il I 


(5.13) 


Upstream, the SPL increased from 0.102 to 0.209, an increase of 6.2 dB. Referring to 
Table 5.6, it is clear that in all cases perturbing one of the design variables so as to 
decrease to sound level downstream produces an increase in sound upstream. Thus, 
it is difficult in this case to simultaneously reduce the radiated noise both upstream 
and downstream. Some physical insight may help to explain this behavior. 

At this point, it may be unclear why changing the blade shape has any effect on 
the radiated noise. Consider subsonic flow over an unloaded fiat plate. Smith [10] 
showed that the magnitude of the outgoing pressure waves may be expressed by an 
equation of the form 

|p| = [ fp(x)}K(M-oo,G,0,uj,<T,x)dx (5.14) 

Jo 

where |p(a:)J is the jump in pressure at some location x on the blade surface, and 
K is a kernel function relating the blade load at x to the outgoing pressure. A 
similar expression may be used to describe the unsteady lift on the blade. Note that 
the kernel function is essentially only dependent on the form of the incoming gust 
and the steady flow conditions upstream of the airfoil. Hence, the kernel function is 
essentially independent of the blade shape. 
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Figure 5.25: Contours of unsteady pressure for cascade of redesigned airfoils due to 
incoming vortical gust at lxBPF. u> = 3.77, a = —158.4°. 
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Next, consider the discretized vector form of Eq. (5.14), i.e., 
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Here the left-hand side of Eq. (5.15) represents the blade row response to the gust, 
consisting of the upstream and downstream going pressure waves and the unsteady- 
lift. The matrix [X] is the discretized form of the kernel function K. The vector T 
represents the discretized values of the pressure jumps along the blade surface. As- 
suming that the duct geometry, upstream flow conditions, and the form of the gust 
are fixed in the design process, the matrix [X], like the kernel function, is essentially 
independent of the blade shape. Therefore, the magnitude of the outgoing pressure 
waves (and the unsteady lift) is only a function of the distribution of the pressure 
jumps on the blade surface. Because changing the blade shape will change the distri- 
bution of the pressure jumps, the magnitudes of the pressure waves will change as a 
function of the airfoil shape. 

The relationship between the outgoing pressure waves and blade loading may also 
help explain why it appears to be difficult to simultaneously reduce the magnitude of 
the upstream and downstream going pressure waves. The response of a cascade to an 
incoming gust generates acoustic radiation in the flow in the form of pressure waves. 
Examination of Figs. 5.17 and 5.23 shows that the magnitude of the unsteady lift is 
nearly the same for the nominal and redesigned airfoil. The acoustic radiation in the 
flow, therefore, should be nearly the same in both cases. Clearly, some of the radiation 
contained in the nominal propagating pressure waves has been transferred from the 
downstream going wave to the upstream going wave. For net noise reduction to be 
achieved, however, some of the radiation must be transferred into acoustic modes that 
do not propagate. In the present two-dimensional analysis, the energy would have 
to be transferred to higher Fourier modes in the circumferential direction, since the 
radial mode shapes are assumed to be uniform. Examination of the magnitudes of 
the pressure waves before and after the redesign shows that some net noise reduction 
has been accomplished. 

The above discussion does not tell the complete acoustic story, however. First, the 
analysis is based an unloaded flat plate theory, while the case in question has steady 
loading and flow turning. Still, the analysis should be applicable in a qualitative sense. 
Second, and perhaps more importantly, all radial effects have been neglected. Chapter 
6 contains a discussion of the three-dimensional acoustic behavior of annular ducts. 
The radial behavior of the pressure modes may be described by a series of orthogonal 
mode shapes, each of which has an associated axial wavenumber that increases with 
the order of the radial mode. Typically, if the pressure waves associated with a 
given Fourier mode propagate, the energy is contained in a small number of these 
radial modes. The higher order radial modes are cut-off. Hence, net noise reduction 
may be achieved by transferring energy from the lower, propagating radial modes to 
the higher, nonpropagating modes. Therefore, a three-dimensional analysis of the 
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Table 5.7: Computational times for present method using 129 x 49 node grid. 


Procedure 

CPU Time (sec) 

Nominal Steady 

837.0 

Nominal Unsteady 

12.8 

Steady Sensitivity Analysis (8 var) 

28.8 

Unsteady Sensitivity Analysis (10 var) 

61.0 

Finite Difference Steady Sensitivity Analysis (8 var) 

13394 

Finite Difference Unsteady Sensitivity Analysis (10 var) 

256 


duct acoustics appears to be important to correctly predicting the radiated noise. 
Furthermore, shape changes along the span of the blade may have a dramatic effect 
on the outgoing pressure. 

5.2.5 Computational Efficiency 

Finally, a note about the computational times for this case. Table 5.7 shows the 
required computational times for the present method and for a finite difference calcu- 
lation. Note that the steady flow calculation times are substantially higher than those 
shown in Table 5.4. In this case, the grid and steady flow equations are solved simul- 
taneously, which increases the cost of the steady flow calculation. The time required 
for the unsteady flow calculation is somewhat higher due to the cost of calculating 
the rotational velocity, v R . Still, the nominal computational time is significantly less 
than a comparable Euler calculation, while retaining all the dominant physics of the 
problem. Comparing the computational times for the sensitivity analysis to the finite 
difference calculation, there is an even more dramatic difference in computational cost 
than the difference shown in Table 5.4. This indicates that the computational savings 
of the present sensitivity analysis greatly increases as the flow model becomes more 
complicated (e.g., the Euler or Navier-Stokes equations). 


5.3 Modern Fan Exit Guide Vane 

The final test case is a three-dimensional fan exit guide vane from a modern high- 
bypass ratio engine. In the previous example, we examined the acoustic response of 
the blade row at a single representative radius of the machine. Furthermore, in the 
previous example we attempted to reduce the downstream acoustic response using all 
of the available design variables. In this example, the objective will be to reduce the 
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Figure 5.27: Contours of magnitude of unsteady pressure for radial station near 
the hub of a modern fan exit guide vane due to incoming vortical gust at 2xBPF. 
iWloo = 0.32, G = 0.46, w = 12.5, a = -288°. 

overall acoustic response using only variables that affect the airfoil shape (i.e., the 
blade-to-blade gap does not change). 

5.3.1 Nominal Analysis 

We will examine the acoustic response and redesign the blade shape at three 
representative radial locations near the hub, midspan, and tip of the vane. Because 
of the solid body rotation in the steady flow due to the upstream row of rotor blades, 
the steady and unsteady flow conditions vary along the span. For brevity, we will not 
examine the steady flow in detail. 

Consider the unsteady flow due to rotor wakes impinging on the EGVs as in 
the previous section. At all radial stations, the lxBPF pressure waves are cutoff. 
Figures 5.27-5.29 show contours of the magnitude of the unsteady pressure due to 
the vortical wakes at 2xBPF. At the radial location near the hub, the inflow Mach 
number, M_oo, is approximately 0.32, the blade-to-blade gap G is about 0.46, the 
reduced frequency, u is 12.5, and the interblade phase angle, a, is —288°. Because 
the angular gap between the blades is constant along the span, the blade-to-blade gap 
increases along the span of the EGV, and the lowest value is at the hub. The Mach 
number in the stator frame varies, due to the relative motion of the rotor blades, from 
approximately 0.32 at the hub to a maximum of about 0.56 near the midspan and 
0.40 near the tip. In a similar fashion, the reduced frequency has a maximum value 
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Figure 5.28: Contours of magnitude of unsteady pressure for radial station near the 
midspan of a modern fan exit guide vane due to incoming vortical gust at 2xBPF. 
oo = 0.56, G = 0.71, u = 7.7, a = -288°. 


near the hub of about 12.5, a minimum value near the midspan of about 7.7, and a 
value of about 10.4 near the tip. 

Note that near the hub the magnitude of the downstream pressure wave is signif- 
icantly higher than the upstream wave. At the midspan, the two waves are nearly 
equal in magnitude, while at the tip the downstream wave again has a larger mag- 
nitude than the upstream wave. Hence, we will primarily attempt to reduce the 
magnitude of the downstream pressure wave. 

5.3,2 Sensitivity Analysis 

Tables 5.8-5.10 show the sensitivity of the steady lift in the cascade direction and 
the outgoing pressure waves to perturbations in ten design variables. Note that the 
steady turning increases along the span so that the steady flow is nearly aligned with 
the axial direction after it passes through the EGVs. For the most part, the steady lift 
sensitivities are very similar from hub to tip, although the actual magnitudes differ. 
In particular, as was the case in the previous exit guide vane analysis, the steady 
loading is strongly influenced by changes in camber, reflex, blade-to-blade gap, and 
stagger. 

The unsteady sensitivities are not as easily characterized. From the “constrained” 
sensitivities, it is clear that thickness and camber play a significant role in the acoustic 
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Figure 5.29: Contours of magnitude of unsteady pressure for radial station near 
the tip of a modern fan exit guide vane due to incoming vortical gust at 2xBPF. 
M_ 0 o = 0.40, G = 0.90, w = 10.4, cr = -288°. 

response of the EGV at all three radial locations. Also, changing the leading edge 
radius and maximum thickness location of the blade results in a moderate change 
in the acoustic response, considering that the leading edge radius has virtually no 
effect on the steady loading. As a result, the “unconstrained” and “constrained” 
sensitivities due to changes in maximum thickness location and leading edge radius 
are nearly the same. 

Note that there does not appear to be a clear pattern to indicate the reason for 
the “trade-off” phenomena between the upstream and downstream pressure waves 
encountered here and in the previous section. For example, changes in camber at the 
hub and midspan require some trade-off between the two pressure waves, while at the 
tip, the magnitude of both pressure waves may be reduced. The lack of consistent 
behavior along the span may further indicate that three-dimensional effects play an 
important role in the overall acoustic response of a blade row, as discussed in the 
previous section. 
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Table 5.8: Change in steady and unsteady flow quantities due to unit perturbations 
in ten design variables near the hub of a modern fan exit guide vane. The nominal 
lift in the cascade direction, L y , is 0.4425, the nominal magnitude of the upstream 
pressure wave, |p up |, is 0.099, and the magnitude of the downstream pressure wave, 
|Pdown|> is 0.350. 




Unconstrained 

Constrained 

Design Variable 


|Pup|' 

IPdowrt | 

IPupl' 

(Pdown | 

Thickness 

0.0253 

0.4339 

0.3562 

0.4292 

0.3937 

Max Thick. Location 

0.0009 

0.0353 

0.0133 

0.0351 

0.0147 

Camber 

1.1135 

0.1606 

1.0434 

-0.6320 

3.9260 

Max Camber Location 

0.2333 

-0.1498 

-0.0747 

-0.0115 

-0.1116 

Leading Edge Radius 

0.0000 

0.0010 

-0.0018 

0.0010 

-0.0018 

Stagger 

-0.3880 

0.0591 

0.2987 

— 

— 

Gap 

0.5000 

-0.0607 

0.3196 

-0.0526 

0.2553 

Reflex 

-1.7127 

0.5068 

0.8007 

— 

— 

Frequency 

— 

0.0341 

-0.0503 

0.0341 

-0.0503 

Interblade Phase Ang. 

— 

0.0609 

-0.5775 

0.0609 

-0.5775 


5.3.3 Redesign of EGV for Reduced Acoustic Response 

Using the sensitivities contained in Tables 5.8-5.10, we now attempt to reduce the 
acoustic response using changes in the airfoil shape. As noted earlier, in this example 
we will only change the airfoil shape — changes in the blade-to-blade gap will not be 
considered. 

At the station near the hub, the leading edge radius was substantially increased, 
along with a small increase in the maximum camber location. The thickness, camber, 
and maximum thickness location were all decreased slightly. Figure 5.30 shows con- 
tours of unsteady pressure for the redesigned airfoil near the hub. In this case, the 
downstream SPL was reduced by 1.0 dB at the expense of a 2.0 dB increase in the up- 
stream SPL. Because the magnitude of the downstream pressure wave is significantly 
larger than the upstream wave, we consider this to be a reasonable trade-off. 

At the station near the midspan, the leading edge radius and maximum camber 
location were also increased, as well as the thickness. The camber and maximum 
thickness location were decreased. Figure 5.31 shows contours of unsteady pressure 
for the redesigned airfoil near the midspan. In this case, net noise reduction was 
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lpl = 0.312 
-1.0 dB SPL 


Figure 5.30: Contours of magnitude of unsteady pressure for radial station near the 
hub of a redesigned modern fan exit guide vane due to incoming vortical gust at 
2xBPF. M_oo = 0.32, G = 0.46, u = 12.5, a = -288°. 



Ipl =0.181 
-0.9 dB SPL 


Figure 5.31: Contours of magnitude of unsteady pressure for radial station near the 
midspan of a redesigned modern fan exit guide vane due to incoming vortical gust at 
2xBPF. M_ oo = 0.56, G = 0.71, cu = 7.7, a = -288°. 
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Table 5.9: Change in steady and unsteady flow quantities due to unit perturbations in 
ten design variables near the midspan of a modern fan exit guide vane. The nominal 
lift in the cascade direction, L y , is 0.6579, the nominal magnitude of the upstream 
pressure wave, |p up |, is 0.209, and the magnitude of the downstream pressure wave, 
|Pdown|) is 0.200. 




Unconstrained 

Constrained 

Design Variable 

Ly 

|Pup|' 

IPdown | 

bupl' 

bdown | 

Thickness 

0.0843 

-0.2852 

-0.1564 

-0.3208 

-0.0647 

Max Thick. Location 

0.0077 

-0.0262 

0.0106 

-0.0294 

0.0190 

Camber 

2.2503 

-0.3106 

0.9183 

-2.3070 

6.5816 

Max Camber Location 

0.3090 

-0.2830 

0.0384 

-0.1242 

-0.5140 

Leading Edge Radius 

0.0000 

-0.0023 

-0.0004 

-0.0023 

-0.0004 

Stagger 

-0.8514 

-0.1605 

0.3151 

— 

— 

Gap 

0.4932 

-0.2381 

0.2108 

-0.2155 

0.1526 

Reflex 

-3.5867 

-0.2646 

0.0623 

— 

— 

Frequency 

— 

-0.1082 

-0.0319 

-0.1082 

-0.0319 

Interblade Phase Ang. 

— 

-0.7778 

-0.0655 

-0.7778 

-0.0655 


achieved, by 0.9 dB downstream and 2.6 dB upstream. 

At the station near the tip, again the leading edge radius and maximum camber 
location were increased. Also, the thickness and camber were slightly increased, 
and the maximum thickness location was decreased. Figure 5.32 shows contours of 
unsteady pressure for the redesigned airfoil near the tip. A small trade-off was made 
here, decreasing the downstream SPL by 1.4 dB while increasing the upstream SPL 
by 0.6 dB. 

Finally, it is clear that the reduction in acoustic response is not as dramatic as 
the reduction achieved in the previous example. One reason for the small reduction is 
that the blade-to-blade gap was held fixed in this example. Although the sensitivities 
shown in Tables 5.8-5.10 indicate that the gap only plays a moderate role in the 
acoustic response, the sensitivities alone do not completely illustrate the role the 
gap plays in the redesign process. Changing the gap also moderates the constraint 
relationship between the camber and reflex so that larger changes in camber may 
be employed, thereby achieving greater noise reduction. Changing the airfoil shape 
alone appears to have only a moderate effect on the radiated noise, at least in the 
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Table 5.10: Change in steady and unsteady flow quantities due to unit perturbations 
in ten design variables near the tip of a modern fan exit guide vane. The nominal 
lift in the cascade direction, L y , is 0.8140, the nominal magnitude of the upstream 
pressure wave, |p up |, is 0.103, and the magnitude of the downstream pressure wave, 
|Pdown|> Is 0.263. 




Unconstrained 

Constrained 

Design Variable 

L'y 

IPupl' 

| Pdown | 

bu P r 

IPdown | 

Thickness 

0.0921 

1.2871 

-1.6826 

1.2239 

-1.8322 

Max Thick. Location 

0.0122 

0.0782 

0.0534 

0.0698 

0.0335 

Camber 

1.9340 

-0.5630 

-0.2809 

-3.2040 

-6.6747 

Max Camber Location 

0.2967 

-0.0227 

-0.0706 

0.1325 

0.3363 

Leading Edge Radius 

-0.0002 

0.0008 

-0.0025 

0.0009 

-0.0021 

Stagger 

-0.7708 

-0.2338 

-0.5205 

— 

— 

Gap 

0.4326 

-0.8044 

1.4336 

-0.7506 

1.5608 

Reflex 

-2.4583 

0.1682 

0.6038 

— 

— 

Frequency 

— 

-0.2758 

0.2425 

-0.2758 

0.2425 

Interblade Phase Ang. 

— 

-0.8420 

0.5372 

-0.8420 

0.5372 


present two-dimensional analysis. 


5.4 Summary 

Although extensive parametric studies have not been performed, the results presented 
in this chapter point to some possible general trends. It appears that the effectiveness 
of the sensitivity analysis increases with reduced frequency. The redesign of the 
aeroacoustic examples resulted in a much larger change in the unsteady aerodynamic 
behavior than the aeroelastic example. This may be due to the difference in frequency 
(an order of magnitude) between the two cases. At high frequencies, the wavelengths 
of the unsteady waves are shorter, and therefore are more dependent on the airfoil 
shape. 

An interesting result from the first aeroacoustic example was that there appears to 
be a trade-off between noise reduction upstream and downstream. In other words, it 
may be difficult to achieve global noise reduction for a particular operating condition. 
The analysis appears to be better suited to cases where there is a noise problem either 
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Ipl = 0.224 
-1 .4 dB SPL 


Figure 5.32: Contours of magnitude of unsteady pressure for radial station near the 
tip of a redesigned modern fan exit guide vane due to incoming vortical gust at 
2xBPF. M-oc = 0.40, G = 0.90, u = 10.4, a - -288°. 

upstream or downstream of the blade row, but not both. An analysis using subsonic 
flat plate theory described a possible reason for this trade-off between the upstream 
and downstream radiated noise, and concluded that three-dimensional effects may 
play a large role in the accurate prediction of the outgoing pressure. 

In all three cases, the steady loading on the redesigned airfoil was moved forward, 
with an increase in the steady pressure gradient on the surface. This points to the 
need for a viscous version of this analysis, so that the increase in steady aerodynamic 
losses can be quantified. 

Finally, the computational efficiency of the analysis indicates that the computa- 
tional savings greatly increases with the complexity of the flow model. Consequently, 
the results presented here would be even more impressive if a Euler or Navier-Stokes 
analysis had been used. 



Chapter 6 

Application to Three-Dimensional 
Problems 


Although the preceding development is well-suited to the design and analysis of 
two-dimensional cascades, it is not entirely obvious how to analyze an actual three- 
dimensional blade row using this method. This chapter provides a framework for 
analyzing three-dimensional aeroacoustic problems using the present two-dimensional 
procedure. Section 6.1 contains a discussion of the general three-dimensional problem 
and the importance of the annular nature of the geometry. In addition, an analyti- 
cal procedure for calculating the natural acoustic modes of an annular duct will be 
developed. Next, Section 6.2 contains a discussion of how the unsteady pressure at 
the far field of a three-dimensional blade row may be calculated. Finally, in Section 
6.3, a general method for calculating the magnitude of the outgoing three-dimensional 
pressure waves in the far field from the computed two-dimensional unsteady potential 
will be described. Also, the application of the sensitivity analysis to the calculation 
procedure described in this chapter will be discussed. 


6.1 Acoustic Modes in an Annular Duct 

The fluid in an axial flow turbomachine is typically modeled as flow through an 
annular duct. A schematic of such a duct is shown in Figure 6.1. The duct has an 
inner (or hub) radius, r#, and an outer (or tip) radius, rj. A key parameter defining 
the duct geometry is the so-called hub-to-tip ratio, th/t?. As this ratio approaches 
unity, the flow becomes nearly two-dimensional. For lower hub-to-tip ratios, however, 
annular effects become more significant. 

The first step in analyzing the three-dimensional acoustic behavior of a blade row 
is to examine the natural acoustic modes in an annular duct such as the one shown in 
Figure 6.1. For a completely general steady flow, this would be a formidable problem 
to solve analytically. Fortunately, the steady flow field is relatively uniform away 
from the blade row (i.e., the far field). Furthermore, in this analysis, all streamtube 
contraction effects are neglected, i.e., the hub and tip radii are constants along the 
axis of the turbomachine. It is also assumed that there is no radial component of the 
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Figure 6.1: Annular duct geometry. 


steady flow. 

For steady flows without swirl, a classic analysis of acoustic modes in annular 
ducts applied to axial flow turbomachines was performed by Tyler and Sofrin [62]. 
The following analysis, to a large extent, will follow their approach. It is assumed 
that the steady flow is entirely axial, with velocity U . Instead of using a reference 
frame fixed in space, however, it will be convenient to solve for the acoustic modes in a 
reference frame attached to the steady flow field, i.e., the relative steady fluid motion 
is zero. In this case, the unsteady pressure field is governed by the three-dimensional 
wave equation in cylindrical coordinates, i.e., 



1 d 2 p 

C 2 ~di 2 


(6.1) 


where 


id_ \_(P_ jp_ 

dr 2 + r dr r 2 dO 2 dx 2 


Here x is the axial coordinate, 8 is the circumferential coordinate, and r is the radial 
coordinate in the reference frame attached to the fluid. Substituting this expression 
into Eq. (6.1) results in 


d 2 p l dp 1 d 2 p d 2 p 1 d 2 p 

dr 2 ^ r dr + r 2 dd 2 dx 2 C 2 dt 2 


(6.2) 


As was assumed earlier, the small disturbance unsteady flow is harmonic with 
known temporal frequency cu, as yet unknown axial wave number a, and known 
circumferential wave number f3 m . In this analysis, the circumferential wave number 



136 CHAPTER 6. APPLICATION TO THREE-DIMENSIONAL PROBLEMS 


will be referred to by an integer m that denotes the circumferential mode number, 


i.e., 


Pm = 


( a + 2mn)9 

©G 


(6.3) 


where 0q is the angular gap between the blades. For now, the radial behavior of 
the unsteady pressure is still unknown. Hence, the behavior of the unsteady pressure 
may be written in the form 


p(r , 9, x, t) = p(r) exp [jm6 + j ax + jut ] 


(6.4) 


where fi{r) is some as yet undetermined function. Substituting this functional form 
into Eq. (6.2) results in 


d 2 p 1 dp m 2 2 u 2 

~dP "I" J~d£ ~ ~P li ~ a ^ + = ° 


If a new variable a* is defined such that 


« =- 2 -a 


then Eq. (6.5) may be written as 




dr 2 




2 -2 
a r 


m 


! ) p = 0 


(6.5) 


( 6 . 6 ) 


(6.7) 


Equation (6.7) is simply Bessel’s Equation. The solution, therefore, is composed of 
Bessel functions of order m, i.e., 


/z(r) = AJ m (a*r) + BY m (a*r) 


( 6 . 8 ) 


The constants A and B are determined by the boundary conditions on the radial 
behavior of the unsteady pressure. Because the hub and tip are solid surfaces, the 
flow tangency boundary condition requires that radial component of velocity on these 
surfaces must be zero. From the radial momentum equation, then, the radial (or 
normal) component of the pressure gradient at the hub and tip surfaces must be zero, 
i.e., 

dp\ 


and 


dr 

dp 

dr 


= 0 


r H 


= 0 


(6,9) 


( 6 , 10 ) 


r T 


Application of the boundary conditions results in a homogeneous system of equations 
for the constants A and B. Cramer’s rule implies that if the constants A and B are 
not both zero, then 


^Kr„) ^(o'rif) 


f'l-’r?) Tf>'rr) 


= 0 


(6.11) 
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Real Amplitude 


Real Amplitude 



Real Amplitude 


Figure 6.2: Typical radial mode shapes pL mn for annular duct. Hub-to-tip ratio, 
r H / r T = 0.5; solid line, n — 0; dashed line, n — 1; dotted line, n — 2. Left, Fourier 
mode m = 1; center, Fourier mode m — 2; right, Fourier mode m~ 8. 


This transcendental equation has an infinite number of roots a* for every integer 

m. For simplicity, the roots will be ordered by magnitude and denoted by the index 

n. Hence, the radial mode shape fi mn (r) has n zeros in the interval rjj < f < 
rj. The mode shapes are also orthogonal with respect to the weighting function r. 
In addition, since Eq. (6.11) only determines the constants A and B to within an 
arbitrary multiplicative factor, the constants may be chosen so that the radial mode 
shapes are orthonormal with respect to the weighting function r, i.e., 


for all l ^ n, and 



J rr T 
r H 


r dr = 0 


fj-mn r dr ~ 1 


(6.12) 

(6.13) 


for all I = n. 

Figure 6.2 shows some typical radial mode shapes fi mn computed using this anal- 
ysis. As noted above, each mode shape has n zeros between the hub and tip radii. 
Now that the radial mode shape is known, the only remaining variable is the axial 
wavenumber a mn . 

The axial wave number a mn corresponding to the radial mode specified by a* mn 
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may be determined from the roots of Eq. (6.11), and Eq. (6.6), so that 


a: 


2 

77171 



*2 

77171 


(6.14) 


In summary, if the steady flow relative to the frame of reference is zero, the behavior 
of unsteady pressure waves in an annular duct may be determined by the temporal 
frequency u>, the circumferential mode number m, and the hub and tip radii, r h and 
r T . 


6.2 Calculation of Far-Field Unsteady Pressure 

We wish to apply the sensitivity analysis procedure developed in this report to the 
complete three-dimensional flowfield. At first glance, it would appear that one pos- 
sible approach would be to simply solve the three-dimensional governing equations 
of the fluid instead of the two-dimensional versions developed here. To obtain the 
computational efficiency of the present method, however, the nominal steady and 
unsteady flow equations must be solved using LU decomposition, i.e., the discretized 
flow equations must be solved directly. Unfortunately, a direct solution of the three- 
dimensional flow equations is not feasible on current computers. Hence, an alternative 
approach to calculating the three-dimensional flowfield needs to be developed. 

Because even iterative fully three-dimensional unsteady aerodynamic analyses of 
cascades have only recently become available for use on modern workstation com- 
puters, most cascade analyses have used a two-dimensional aerodynamic model. The 
blades are analyzed in what is referred to as “strip-theory.” In this approach, the un- 
steady aerodynamic loads at a number of radial locations along the blade surface are 
computed using a two-dimensional aerodynamic model. In the aeroelastic problem, 
these aerodynamic loads are then integrated along the span to obtain modal forces. 
In the aeroacoustic problem, the unsteady pressure at the far field p(r, x ) may be 
computed from these iaerodynamic loads using a Green’s function approach, i.e., 

p(r,x) = J^lp(r 0 ,x 0 )}G{r,x]r 0 ,Xo)dT (6.15) 

where |p(r<j, x 0 )} is the jump in pressure at the point ( r 0 ,x 0 ) on the blade surface, 
and G is the Green’s function derived for a specific annular geometry and steady flow 
condition [1, 58]. The Green’s function may be considered an “influence coefficient” 
relating the loads on the surface to the pressure in the far field. The Green’s function 
approach relies on a relatively idealized model, however. Typically, the model does 
not include blade thickness or camber, and may only approximate the actual steady 
loading on the blade. 

The linearized potential analysis described in this report improves on the Green’s 
function approach by solving the governing equations of the steady and unsteady 
flow using modern computational fluid dynamic techniques. Arbitrary blade shapes 
may be analyzed, the actual steady blade loading is accounted for, and the far-field 
pressure is computed with the rest of the unsteady flow solution. 
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Although the strip theory approach for modeling quasi-three-dimensional effects is 
a simple and straightforward method, it should be noted that in some cases the results 
computed in this way may be misleading. In particular, Hall and Lorence [35] found 
that real three-dimensional effects may be extremely important in the aeroelastic 
behavior of fans and compressors with low hub-to-tip ratios. 

With this in mind, the next step is to calculate the magnitudes of the actual 
three-dimensional outgoing pressure waves in the duct using the two-dimensional 
analysis developed in this report. Specifically, the two-dimensional unsteady potential 
calculated at the far field will be matched with a three-dimensional analysis of the 
surrounding duct. 


6.3 Calculation of the Outgoing Pressure 


As described earlier in this chapter, Tyler and Sofrin developed an analytical descrip- 
tion of the small disturbance behavior of acoustic modes in annular ducts. However, 
their analysis is limited to flows without mean flow swirl. For more general steady 
flows, a more complicated model is required. For example, Kerrebrock [59] has de- 
veloped an analytical model for certain cases of swirling flows. General analytical 
solutions for flows in annular ducts, however, are not available. If the mean flow field 
is axisymmetric, however, the shape of the eigenmodes in the circumferential direc- 
tion will be Fourier modes. Hall, Lorence, and Clark [51] developed a procedure to 
take advantage of this behavior in their formulation of fully three-dimensional non- 
reflecting boundary conditions for the linearized Euler equations. Their method will 
be adapted here to compute the magnitude of the actual three-dimensional outgoing 
pressure waves using the unsteady potential computed from the analysis developed 
in this report. 

The modeling of the three-dimensional flowfield begins with the full unsteady 
Euler equations in cylindrical coordinates expressed in the rotating frame of reference, 
which may be written as 


dll dF ldG 1 dr H - 

~Tu "n I n7T a S — 0 

Ot Ox r do r or 


(6,16) 


Here x, 6 , and r are the coordinates in the axial, circumferential, and radial directions. 
Also, U is the vector of conservation flow variables, and F, G, and H are the flux 
vectors, and S is a source term due to rotation. These terms are given by 
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pv 
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pw 


0 

puw 


0 

pvw 

s = 

—pw(i>/r — 20) 

pw 2 + p 


p(v — Or) 2 jr + pjr 

pwl 


0 


where p is the density, p is the pressure, e is the interna] energy, and I is the rothalpy. 
Here we have assumed that the coordinate system is rotating about the x-axis with 
rotational speed O. The pressure, p, and the rothalpy, / are given by 


P = (7 - 1) 


e — \-p{u 2 + v 2 + w 2 ) 4- \pVt 2 r 2 

A 


and 


I = 


e + p 


= ^lf + 5 (i! + i3 + “ 2) -5 nv 


where r is the distance from the x-axis (r = \Jy 2 -j- z 2 ). Note that it, v, and w are now 
the flow velocities in the axial, circumferential, and radial directions, respectively. 

It is assumed that the mean flow, U, is axisymmetric and uniform in the axial 
direction. The small disturbance behavior in the far field is then governed by 


<9u dE du 1 dG du 1 d ( rd H \ dS 
~dt + dXJdH + r&U~d9 + rlfr \dU U ) ~ 3U U = ° 


(6.17) 


where u is the perturbation flow and dE/dU, dGjdD , dH/d\J, and dS/dXJ are 
Jacobians based on the mean flowfield, U. 

As discussed in Chapter 2, it is assumed that the shape of the eigenmodes in 
the circumferential direction are Fourier modes. In addition, it is assumed that the 
frequency of the unsteady excitation is u and the interblade phase angle is cr. Using 
these assumptions, the unsteady conservation variables may be expressed as the series 

OO OO 

u (x,0,r,i)= £ £ w mn u mn (r) e j ^ 6 ^ (6.18) 


where (3 m = (cr + 2i rm) is the circumferential wave number of the mth Fourier mode, 
Oq is the angular gap between any two adjacent airfoils, w mn is a coefficient which 
indicates how much of a given mode u mn is contained in the solution u. 

Substituting Eq. (6.18) into Eq. (6.17), and noting that each term in the series 
must vanish gives 


[!] ^mn + ja v 


' dF ' 

] ftm 

dG' 

d\J 

Urnn r\ 

rO a 

dv _ 


u r 


1 d 

A r 

r dr V 


<9H 


dlJ 




dS 

dV 


Umn — 0 


(6.19) 
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Equation (6.19) is an eigenvalue problem for the mnth eigenvalue a mn and associated 
radial eigenvector u Tnn (r). Because the hub and tip casings are solid surfaces, this 
eigenvalue problem is solved subject to the condition that there be no radial flow on 
these surfaces. Since no closed form solution exists for Eq. (6.19) (except for a few 
special cases), the eigenmodes and eigenvalues must be found numerically. This may 
be accomplished by discretizing Eq. (6.19) on a one-dimensional computational grid 
which extends from the hub to the tip of the duct at the far-field boundary. The 
derivatives in Eq. (6.19) may be approximated using finite difference operators. For 
each Fourier mode m, the resulting eigenvalue problem will be of the form 

[M m ]u mn = a mn [N m ]u TOn (6.20) 

where [M m ] and [N m ] are sparse, complex, non-Hermitian matrices of size 5 K x 5 K 
where K is the number of grid points in the radial direction. 

To determine the direction the eigenmodes travel, we examine the wave numbers 
of the eigenmodes. As described in the far-field boundary condition section of Chapter 
2, waves that have complex wave numbers decay as they propagate and hence do not 
contribute to the radiated noise. We only wish to examine the waves which propagate 
unattenuated (i.e., the wave number is purely real). Recall from Chapter 2 that the 
direction of propagation is determined by the group velocity, given by Eq. (2.96). 
Because the propagating pressure waves have distinct eigenvalues, we may calculate 
da/du using the expression [60] 

doi mn 

duj 

where v^ n and u mn are the left and right eigenvectors of the mode. Hence, once the 
eigenvalues and left and right eigenvectors of the individual modes are known, the 
group velocity may be determined, which specifies the direction of propagation. 

The next step is to determine the far-field values of the conservation variables from 
the two-dimensional analysis. Because the steady flow is subsonic, there is a simple 
linear relationship between the primitive and conservation variables, i.e., u = £[u p ], 
where u p is the vector of primitive variables and £ is a linear operator. Hence, for 
simplicity, we write the primitive variables in terms of the unsteady potential, <f>, so 
that 

p - -g${jw<f> + V$ • V<£) 

u - dfyjdx 

Up = < v = d(j)/dy ‘ (6.22) 

w = 0 

p = — R (ju>(f> + V$ • V<^) 

Note that the radial component of the velocity, w, is zero because the potential has 
been calculated using a two-dimensional aerodynamic model. 


,T [ 9[M m ] MNmIl u 

'run dui mn du r 




(6.21) 
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Next, the unsteady conservation variables are Fourier transformed at each radial 
station at the far field. The mth Fourier coefficient is given by 

r&G 


U r 


&G JO 


uexp [-j/3 m 0/Q G ] d9 


(6.23) 


Equation (6.23) may be approximated numerically using the trapezoidal rule. 

For each Fourier mode, in general the vector u m is composed of both incoming 
and outgoing eigenvectors, u mn . We wish to modify the solution (now denoted u^ d ) 
so that only the contribution due to the outgoing propagating modes are retained. To 
accomplish this, we discretize Eq. (6.19) to obtain the eigenmode description of the 
mth Fourier mode. Having computed these mode shapes, we expand u^ d as a sum of 
eigenvectors as in Eq. (6.18) Using the biorthogonality condition, the coefficient w mn 
may be expressed as 


'UJmn — 


JNmlu: 
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old 
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(6.24) 


Hence, to eliminate all but the outgoing propagating modes from the solution, we let 


u_ 


= £ u 


m n 

77171 rj\ 


u: 


jNr, 

[N m ]u r 


old 


= [Z m ]u 


old 


(6.25) 


In this form, Eq. (6.25) acts like a linear filter, removing all but the propagating 
outgoing modes. This process is repeated for each Fourier mode. The modified 
Fourier coefficients u£f w are then summed using Eq. (6.18) to form the new solution 
u at the far-field boundary. The new solution may then be used to calculate the 
sound pressure level (SPL), sound power level (PWL), or other measure of the sound 
output from the blade row. The process may then be repeated for other multiples of 
the blade passing frequency (BPF) to determine the complete tonal noise signature 
of the blade row. 

Equation (6.25) is an important result. Recall that the computation of the eigen- 
analysis required only the duct geometry and the steady flow at the far field. Hence, 
the matrix [Z m ] is not dependent on the airfoil shape if the steady flow angle at the 
far field does not change. 

The calculation of the perturbation of the radiated tonal noise, then, is actually 
quite straightforward. From Eq. (6.25), the perturbed solution at the far field may 
be written as 

uT = [Z™)u^ d - [Z m ]'< (6.26) 

where the primes refer to the perturbed quantities. If the duct geometry and steady 
flow at the far field do not change (a very reasonable assumption), then [Z m ]' is 
zero, and there is simply a linear relationship between the perturbation of the solu- 
tion calculated by the two-dimensional analysis and the perturbed three-dimensional 
radiated noise. 

In summary, because of the simplicity of the model, this method is an efficient 
approach to more realistic three-dimensional modeling of the aeroacoustic behavior 
of a blade row. Furthermore, the method is well-suited to the sensitivity analysis 
developed in this report, and is consistent with three-dimensional analysis techniques 
currently in use. 



Chapter 7 

Conclusions and Future Work 


This chapter summarizes the important features of the present method, and includes 
the important conclusions from the computed results. In addition, some suggestions 
for future work in this area are offered. 


7.1 Conclusions 

This report has presented a new method for calculating the aeroacoustic and aeroelas- 
tic response of a cascade to small changes in the airfoil or cascade geometry. Previous 
methods have calculated the steady and unsteady aerodynamic behavior of a cascade 
based on the airfoil shape and flow conditions. The present method addresses how 
the steady and unsteady flowfields are dependent on the airfoil shape and cascade ge- 
ometry. The information from this sensitivity analysis can be used to suggest design 
changes for improved aeroacoustic and/or aeroelastic performance. In this report, 
the procedure has been applied to the linearized potential equation. The method is 
general in nature, however, and may be applied to other flow models. 

The present method begins with the assumption that the nonlinear unsteady flow 
may be split into a nonlinear steady (or mean) flow, and a small disturbance un- 
steady flow that is harmonic in time. The steady flow is modeled by the steady form 
of the full potential equation, which is discretized using a variational finite element 
technique. The discretized steady equations are solved using Newton iteration with 
LU decomposition at each iteration. Once the steady flow has been computed, the 
unsteady flow is governed by a set of variable coefficient differential equations (the 
discretized form of the linearized potential equation), where the coefficients are de- 
pendent on the steady flow. The linearized potential equation is also discretized by 
a variational finite element procedure. For flutter problems, the variational principle 
incorporates a deforming grid to improve the accuracy of the unsteady solution. For 
forced response problems, the variational principle includes the effect of incident vor- 
tical gusts using rapid distortion theory. The discretized small disturbance equations 
are solved using a single LU decomposition. 

The sensitivities of the steady and unsteady flow to changes in design variables 
are computed in a straightforward fashion by expanding the steady and unsteady flow 
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equations in perturbation series and collecting terms of equal order. The resulting 
equations for the sensitivities may be solved very efficiently by performing one forward 
and one back substitution for each design variable using the LU factors obtained from 
the nominal solutions. 

To demonstrate the effectiveness of the procedure, representative examples of 
aeroelastic and aeroacoustic problems in turbomachinery were presented. In the first 
example, a cascade of compressor blades that was aeroelastically unstable in torsion 
was redesigned to be aeroelastically stable. In the second example, a fan exit guide 
vane was redesigned to significantly reduce the downstream acoustic response. In the 
final example, a fan exit guide vane was redesigned to reduce the overall acoustic 
response using only changes in the blade shape. In all three cases, the redesign 
was accomplished while keeping the steady work done by the cascade and the airfoil 
incidence angle unchanged. 

Analysis of the computed results indicates that the sensitivity analysis gives excel- 
lent predictions of the the actual redesigned airfoil response. In addition, it appears 
that the effectiveness of the modifying the blade shape to improve unsteady aero- 
dynamic performance increases with frequency, i.e., high frequency problems appear 
to be more heavily dependent on the details of the airfoil shape. For aeroacoustic 
problems, there appears to be a trade-off between upstream and downstream acoustic 
radiation, indicating that it may be difficult to simultaneously reduce the upstream 
and downstream radiated noise, at least in the present two-dimensional analysis. Fi- 
nally, the sensitivity analysis was found to be computationally very efficient, with the 
sensitivity analysis requiring a fraction of the time necessary for the nominal calcu- 
lation, and orders of magnitude less time than a finite difference calculation would 
require. Furthermore, the computational efficiency dramatically increases with the 
size of the problem, so the savings will be even more dramatic for flow models with 
increased sophistication. 

7.2 Future Work 

Although the results presented here are extremely encouraging, more work needs to 
be done to fully realize the potential of this method as a design tool. This section 
will discuss some possible avenues for future research. 

7.2.1 Other Flow Models 

Throughout this report, it was emphasized that the sensitivity analysis procedure 
presented here is general in nature and may be applied to a wide variety of flow 
models and discretization schemes. Hence, the first logical extension of the present 
method is to apply it to a more sophisticated flow model. 

Linearized Euler analyses of unsteady flow problems have become more common 
in recent years [30, 31, 35]. Unlike the linearized potential formulation, the Euler 
equations allow transonic problems including shocks to be modeled accurately by 
including the effect of vorticity and entropy generation. One advantage of this ap- 
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proach is that a sensitivity analysis using the linearized Euler equations could be used 
to predict the movement of the shock location due to changes in the airfoil shape. 

A more powerful application would be to apply the sensitivity analysis proce- 
dure to the linearized Navier-Stokes equations. Recently, Clark and Hall [34] used 
the linearized Navier-Stokes equations to examine stall flutter in cascades. Although 
they used an iterative approach to calculate the steady and unsteady solutions, their 
procedure could be used to generate the matrix entries for a direct solution method. 
There are several advantages in using a Navier- Stokes-based approach. For example, 
the steady losses due to blade shape changes can be quantified. This would allow the 
designer to examine the actual trade-off between steady efficiency and unsteady per- 
formance. Another example is that the sensitivity analysis could be used to examine 
stall flutter problems. Specifically, the analysis could suggest design changes to help 
prevent flow separation, or at least blade design changes which suppress stall flutter. 
Although Navier-Stokes-based analyses are relatively immature, a sensitivity analysis 
based on the Navier-Stokes equations would be a powerful design tool. 

7.2.2 Multidisciplinary Optimization 

As was demonstrated in Chapter 5, the aeroacoustic and aeroelastic design strategy 
currently implemented involves simply using the results of the sensitivity analysis di- 
rectly. A more sophisticated design strategy would be to consider finding the optimal 
design of an airfoil for an operating condition. For example, a designer may wish 
to minimize the radiated noise subject to the constraint that the desired turning is 
achieved. Additional constraints may be necessary. For example, a constraint on the 
surface pressure gradient may be imposed so that the flow does not separate. The 
resulting optimal solution may be determined through the use of a nonlinear con- 
strained minimization procedure [63]. The problem is nonlinear because the acoustic 
behavior of the cascade is a nonlinear function of the blade shape. Efficient opti- 
mization procedures require accurate estimates of the sensitivity of the cost function 
being optimized to changes in the design variables. The previously described sensitiv- 
ity analysis is well-suited to such a procedure because the sensitivities are calculated 
accurately and efficiently. 

Another advantage of an optimal design procedure is that it may be multidisci- 
plinary. The sensitivity analysis procedure developed in this report only addresses 
the aerodynamic behavior of the aeroelastic or aeroacoustic problem under consider- 
ation. Similar analyses could be performed, for example, on the equations governing 
the blade structure. In the optimization procedure, both the structural and aero- 
dynamic sensitivities could be combined to create a multidisciplinary optimal design 
procedure. 

7.2.3 Multiple Blade Rows 

In this report, the blade row being analyzed is assumed to be isolated in an infinitely 
long duct. In reality, of course, these blade rows are not isolated, and in fact the blade 
rows in a turbomachine are spaced together quite closely. As a result, it is necessary 
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to model the blade row interaction to obtain a more complete model of the flowfield 
within a turbomachine. Unfortunately, blade row interaction is a very complicated 
problem to mode! because of the relative motion between blade rows. This relative 
motion produces a shifting of the frequencies of the acoustic, vortical, and entropic 
waves as they propagate from one blade row to the next. Furthermore, the blade 
rows reflect these waves at a number of shifted frequencies and scattered interblade 
phase angles, violating the assumption that the solution is harmonic in time with a 
fixed interblade phase angle. 

Recently, Hall and Silkowski [64] and Hanson [65] have developed a technique 
to analyze blade row interaction by analyzing each blade row independently, and 
computing the reflection and transmission coefficients for each pressure and vorticity 
mode in the duct. A set of linear equations may then be constructed that couples 
together these coefficients to determine the interaction. Although these methods so 
far have used only flat-plate unsteady aerodynamics, in principle any linearized so- 
lution procedure could be used to model the unsteady flowfield. Hence, the present 
sensitivity analysis could be applied to this procedure. Changes to the airfoil shape 
or cascade geometry will result in perturbations in the reflection and transmission co- 
efficients associated with the blade row. Also, perturbations in the incoming pressure 
and vorticity waves to a blade row will also affect its unsteady response. Once these 
sensitivities are calculated, the coupling equations may be expanded in a perturbation 
series. The first-order terms will determine the change in the blade row interaction 
due to a change in the blade shape. 

7.2.4 Three-Dimensional Problems 

In Chapter 6, some aspects of the three-dimensionality of unsteady flows in tur- 
bomachines were examined. A framework was then presented of how the present 
two-dimensional analysis could be used to analyze three-dimensional problems. Al- 
though the approach described in Chapter 6 is clearly better than simply ignoring 
three-dimensional effects, it does not model fully the annular nature of the blade 
row. To do this, a fully three-dimensional model is required. Unfortunately, most 
three-dimensional unsteady aerodynamic models are extremely computationally ex- 
pensive, due to the large number of computational nodes required. Potential methods 
are generally not suitable for three-dimensional analyses because they are unable to 
model the mean flow swirl produced by the rotating blades. As a result, the Euler 
equations model the minimum amount of physics necessary to realistically examine 
the three-dimensional unsteady aerodynamic problem. Three-dimensional linearized 
analyses are relatively new, the first implementation being the Euler method of Hall 
and Lorence [35]. 

Unfortunately, even if a linearized solver has been obtained, extending the pre- 
vious two-dimensional approach to three dimensions is not a straightforward task. 
In particular, LU factorization requires far too much computer time and computer 
storage to be feasible. To overcome this problem, Hall and Lorence used an iterative 
method to compute the steady and unsteady flowfields. This approach, however, is 
not useful for the sensitivity analysis because the computation of the sensitivities 
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would take as much time as the nominal calculation. 

One possible approach to the three-dimensional sensitivity analysis would be to use 
a preconditioned conjugate gradient (CG) method, such as GMRES [66]. Burgreen 
and Baysal [67] have used such an approach for steady aerodynamic optimization 
problems. In addition, to reduce the computational expense, they used the domain 
decomposition scheme developed by Eleshaky and Baysal [68], where the computa- 
tional domain is divided into subdomains, each of which contains internal cells and 
boundary/interface cells. In this method, the subdomains are solved independently 
after the influences of the boundary cells is calculated. Unfortunately, the GMRES 
algorithm would still have to be applied for each right-hand side, which significantly 
reduces the efficiency of the sensitivity analysis. 

A solution method developed for solution of a system of linear equations with 
multiple right-hand sides is required. Consider the following approach. Instead of 
insisting on exactly factoring the large sparse matrices in the steady and unsteady 
solution procedure, an approximate but accurate factorization based on the Lanczos 
algorithm may be sufficient. The Lanczos algorithm is a recursive algorithm which 
factors a matrix into a symmetric tridiagonal matrix and an orthonormal matrix con- 
taining the so-called Lanczos vectors at each iteration [69]. In principle, the recursive 
algorithm exactly factors the matrix after n iterations where n is the number of un- 
knowns. However, unlike LU decomposition, the recursion may be stopped after only 
m iterations where m <C n to obtain an approximate factorization of the matrix. 
In this case, only m Lanczos vectors are retained, and the solution to the system of 
equations may be obtained for considerably less computational effort than solving the 
system exactly. Furthermore, the Lanczos vectors only depend on the matrix, not 
the right-hand side, so once the Lanczos vectors have been computed, the solution for 
multiple right-hand sides may be obtained for little additional computational effort. 

The proposed algorithm is similar in many respects to the popular GMRES algo- 
rithm. In fact, the conjugate gradient algorithm may be derived from the Lanczos 
algorithm. One of the disadvantages of the conjugate gradient method, as noted 
earlier, is that it is a “one-shot” method. For each right hand side, the GMRES 
algorithm must be restarted. The Lanczos approach, on the other hand, retains the 
approximately factored matrix for reuse. This approach has the potential to extend 
the encouraging results presented in this report to three dimensions. 
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Sensitivity of the NACA Modified 
Four-Digit Airfoil Definition 


As an example of how airfoil shapes and their perturbations are calculated, the NACA 
modified four-digit series of airfoils will be examined [61]. The NACA four-digit 
definition allows the specification of five variables to determine the airfoil shape: the 
magnitude of maximum thickness, m t , the magnitude of the maximum camber, m c , 
the chordwise location of the maximum thickness, l u the chordwise location of the 
maximum camber, 4, and the leading edge radius, r t . 

The thickness distribution is defined by the following two equations: 

± yt — a 0 y/x + a x x + a 2 x 2 + a 3 x 3 (B.l) 

ahead of the location of maximum thickness, and 

± y t = d 0 T cfi(l - x) + d 2 { 1 - x ) 2 + d 3 (l - a;) 3 (B.2) 


aft of the maximum thickness location, where lengths have been nondimensionalized 
by the blade chord, c. The positive sign of y t corresponds to the upper surface of the 
airfoil, while the negative sign refers to the lower surface. 

The four coefficients d 0 , d+, d 2 , and d 3 are determined from four boundary condi- 
tions. The first is that the thickness at the maximum thickness location is one-half 
the magnitude of the maximum thickness. Second, the thickness at the trailing edge 
is zero. Third, the slope of the thickness distribution at the maximum thickness lo- 
cation is zero. Fourth, the slope of the thickness distribution at the trailing edge is 
specified, so that the airfoil has a finite wedge angle at the trailing edge. The specified 
slope may be determined from the following curve fit: 


dy t 

dx 


= -mt (0.775 + 2.51667x - 13.625x 2 + 35.8333x 3 - 12.5x 4 ) (B.3) 

TE 


where TE refers to the trailing edge of the airfoil. 

The four coefficients a 0 , «i, a 2 , and a 3 are also determined from four boundary 
conditions. First, as before, the magnitude of the thickness at the maximum thickness 
location is specified. Second, the slope of the thickness distribution at the maximum 
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thickness location is zero. The third condition is that the curvature (or second deriva- 
tive) of the thickness distribution at the maximum thickness location is the same as 
that computed from the aft coefficients. Fourth, the leading edge radius is speci- 
fied, i.e., the thickness at the leading edge radius location is the leading edge radius 
{yt = x )• According to the definition, the leading edge radius is 

r t = 1-1019 (B.4) 

where I is the first digit following the dash in the designation and the value of I does 
not exceed 8. 

For each set of coefficients, a linear matrix equation may be constructed. The 
solution of these two equations results in the coefficients. Functionally, it is clear that 

d = d (m t ,£ t ) (B.5) 

and 

a = a (m t ,e t ,r t ) (B.6) 

where d and a are vectors containing the defining coefficients. 

The specification of the mean line is considerably more straightforward. The mean 
line is defined by two equations: 

2/m = ( 2 ^ x - * 2 ) ( B - 7 ) 

ahead of the maximum camber location, and 

7Ti l 

y m = (1 ~ ° £ - 2 [(1 - 24) + 2i e x - x 2 } (B.8) 

aft of the maximum camber location. Note that the shape of the mean line is simply 
two parabolas connected together with (7 1 continuity. 

The actual airfoil definition is the superposition of the mean line and the thickness 
distribution, i.e., 

2 /air = y m ±yt (B.9) 

where the surface is determined by the sign of y t . In practice, the airfoil definition 
points are calculated at some finite number of locations, N. This set of ( x , y) locations 
is then rotated through the stagger angle, 0, about a specified rotation point (x r ,y r ), 
so that 

Xi = (XAIRi - x r) COS 0 - (j/AIRi “ 2/r) sin 0 

1,2,..., AT (B.10) 

Y % = ( 2 ; air* - x r ) sin 0 + (y A iR,- - 2/r) cos 0 

where Xi and Y{ are the airfoil defining shape points discussed in Chapter 3. The arc 
length, Si, is then computed at each point, and the X and Y locations are splined 
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as a function of S. The X and Y locations may be used to form the airfoil shape 
portion of the vector Z in the grid and steady flow equations, Eqs. (3.8) and (3.9). 

Now that the airfoil definition and spline generation procedure has been explained, 
the perturbation of the airfoil defining points, Z', may be computed. Consider the 
effect of a small change in the magnitude of the airfoil camber, m' c . The perturbation 
of the unrotated points defining the airfoil shape is 

2/air = yr (% tc x ~ x ) (B.ll) 

ahead of the maximum camber location, and 


2/air — 


( 1 - 4) 2 


[(1 - 2 Q + 2 4* - I 2 


(B. 12) 


aft of the maximum camber location. Because the ^-locations of the defining points 
do not change, :Ea IR is zero. 

In general, then, the perturbation of the y-locations of the defining points may be 
written as: 


, _ dyr 

2 / AIR — 


dm, 


-m' + 


dy r 

dt 




(B.13) 


This perturbation at each point is then rotated through the stagger angle, so that 


X[ - (xaiR; - x ' T ) cos 0 


(2/air; — y' r ) sin 0 


-(zAiRi - x r )0' sin0 + 


(j/AIR i - 2 /r) 0 ' COS 0 


(B.14) 



(zair; ~ x 'r) sin 0 


+ ( 2 /AIR i ~ y'r) COS 0 


+(x A ir i - x r )0' COS 0 


( 2 /AIR i - 2 /r) 0 ' sin 0 


(B.15) 


where 0' is the perturbation in the stagger angle, and x' T and y' T are the perturbations 
in the x- and y-location of the rotation point, which are included here for complete- 
ness. The perturbation in the arc length at each point, S', may also be computed. 
The perturbations are then cubic splined as a function of arc length. 

In summary, then, this appendix has shown how the vector of defining points 
that form part of the vector Z could be calculated for an actual airfoil shape. The 
perturbation in these points, Z', due to small changes in airfoil defining parameters 
may be calculated by perturbing the analytical expressions that define the airfoil. 
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